Water-saving control method and system for agricultural irrigation depth based on multi-source data fusion

Abstract

The application provides an agricultural irrigation depth water-saving control method and system based on multi-source data fusion, which collects multi-dimensional original monitoring data, generates a monitoring data field, extracts physiological rhythm features and environmental fluctuation features under different time dimensions through feature extraction, obtains a crop water sensitivity feature set through feature correlation analysis screening, inputs a crop water metabolism model, generates a water supply and demand balance coefficient, constructs an irrigation decision optimization model to obtain an irrigation parameter set, drives an irrigation execution mechanism to execute a partition variable irrigation operation, collects crop transpiration rate change data after irrigation, corrects the water supply and demand balance coefficient, and generates a final irrigation control instruction to drive the irrigation execution mechanism to execute the partition variable irrigation operation. The application can continuously optimize the irrigation control instruction, ensure the dynamic matching of the irrigation decision and the actual water demand of crops, and effectively improve the water-saving effect and water use efficiency of agricultural irrigation.

Description

A Deep Water-Saving and Water-Controlling Method and System for Agricultural Irrigation Based on Multi-Source Data Fusion Technical Field

[0001] This invention relates to the field of data processing, and more specifically, to a method and system for deep water-saving and water-controlling in agricultural irrigation based on multi-source data fusion. Background Technology

[0002] With the increasing prominence of global water scarcity, the optimization and upgrading of water-saving and water-control technologies in agricultural irrigation, as a major water user, has become crucial for sustainable agricultural development. Deep water-saving and water-control technologies in agricultural irrigation aim to maximize water use efficiency and reduce irrigation water consumption while meeting the water requirements of crop growth through precise regulation of the irrigation process. Current technologies lack systematic integration of crop physiological responses and environmental factors during the data acquisition phase, resulting in insufficient data correlation mining. The feature processing stage struggles to accurately capture the dynamic coupling relationship between crop root water uptake and soil moisture conduction, leading to insufficient matching between irrigation parameters and actual crop water requirements, thus affecting the water-saving effect and water use efficiency improvement of agricultural irrigation. Summary of the Invention

[0003] In view of this, the present invention provides a method and system for deep water-saving and water-controlling in agricultural irrigation based on multi-source data fusion.

[0004] According to one aspect of the present invention, a method for deep water-saving and water-controlling agricultural irrigation based on multi-source data fusion is provided. The method includes: collecting multi-dimensional raw monitoring data of an agricultural production area to generate a monitoring data field with spatiotemporal continuity, wherein the monitoring data field includes a crop physiological response time-series data layer and an environmental factor dynamic data layer; extracting features from the monitoring data field to extract physiological rhythm features and environmental fluctuation features under different time dimensions, and obtaining a crop water-sensitive feature set through feature correlation analysis, wherein the crop water-sensitive feature set includes key feature indicators characterizing changes in crop water status; inputting the crop water-sensitive feature set into a crop water metabolism model to infer the coupling relationship between crop root water absorption rate and soil water conductivity, thereby generating a model characterizing crop water supply and demand. A dynamic equilibrium water supply and demand balance coefficient is required. Based on the water supply and demand balance coefficient and real-time soil moisture data, an irrigation decision optimization model is constructed. The irrigation decision optimization model includes an irrigation total quantity optimization model that optimizes the total irrigation amount with the goal of maximizing water saving rate and a growth stage irrigation quota allocation model that allocates irrigation quotas for each growth stage with the goal of maximizing water use efficiency. An irrigation parameter set is obtained through the irrigation total quantity optimization model and the growth stage irrigation quota allocation model. The irrigation execution agency is driven to perform zonal variable irrigation operations according to the irrigation parameter set. Data on crop transpiration rate changes after irrigation are collected. The water supply and demand balance coefficient is corrected based on the crop transpiration rate change data, and a final irrigation control command is generated to drive the irrigation execution agency to perform zonal variable irrigation operations.

[0005] According to another aspect of the present invention, a computer system is provided, comprising: a processor; and a memory, wherein the memory stores a computer program that, when executed by the processor, causes the processor to perform the method as described above.

[0006] This invention provides a method for deep water-saving and water-controlling agricultural irrigation based on multi-source data fusion. It generates a spatiotemporally continuous monitoring data field by collecting multi-dimensional raw monitoring data from agricultural production areas under a preset spatiotemporal sampling grid, achieving structured integration of crop physiological responses and environmental factor data. By extracting physiological rhythm features and environmental fluctuation features at different time dimensions from the monitoring data field and filtering them to obtain a crop water-sensitive feature set, it can accurately focus on key feature indicators strongly correlated with changes in crop water status, improving the targeting and effectiveness of feature extraction and overcoming the feature redundancy problem caused by experience dependence in traditional feature extraction processes. The crop water-sensitive feature set is input into a crop water metabolism model to infer the coupling relationship between crop root water absorption rate and soil water conductivity. The generated water supply and demand balance coefficient can dynamically reflect the real-time balance state of crop water supply and demand, overcoming the limitation of traditional static parameter analysis in failing to capture the dynamic interaction of the crop-soil water system. Based on the water supply and demand balance coefficient and real-time soil moisture data, an irrigation decision optimization model is constructed. Through the synergistic effect of an irrigation total quantity optimization model and a growth stage irrigation quota allocation model, it achieves irrigation parameter optimization with the goals of maximizing water saving rate and water use efficiency. By driving the irrigation actuators based on the set of irrigation parameters and correcting the water supply and demand balance coefficient based on the changes in crop transpiration rate after irrigation, irrigation control instructions can be continuously optimized, ensuring dynamic matching between irrigation decisions and actual crop water needs, and effectively improving the water-saving effect and water use efficiency of agricultural irrigation. Attached Figure Description

[0007] Figure 1 is a schematic diagram of the application scenario provided by the embodiment of the present invention.

[0008] Figure 2 is a schematic diagram of the implementation process of a deep water-saving and water-control method for agricultural irrigation based on multi-source data fusion provided by an embodiment of the present invention.

[0009] Figure 3 is a schematic diagram of the hardware entity of a computer system provided in an embodiment of the present invention. Detailed Implementation

[0010] The water-saving and water-control method for agricultural irrigation based on multi-source data fusion provided in this embodiment of the invention can be applied to the application environment shown in Figure 1. The monitoring device 102 communicates with the computer system 104 via a network. The data storage system can store the data that the computer system 104 needs to process. The computer system 104 can obtain monitoring data fields from the local storage of the monitoring device 102, the data storage system, or cloud storage. The monitoring device 102 can be, but is not limited to, a time-domain reflectometer, a laser displacement sensor, a humidity sensor, etc. The computer system 104 can be implemented using a standalone server or a server cluster consisting of multiple servers.

[0011] The water-saving and water-control method for agricultural irrigation based on multi-source data fusion provided in this invention is applied to a computer system, and specifically includes the following steps:

[0012] Step S100: Collect multi-dimensional raw monitoring data of agricultural production areas to generate a monitoring data field with spatiotemporal continuity. The monitoring data field includes a crop physiological response time-series data layer and an environmental factor dynamic data layer.

[0013] Specifically, data can be collected under a preset spatiotemporal sampling grid. This grid is pre-defined and used to divide the agricultural production area spatially and temporally to facilitate regular data collection. Spatially, the agricultural production area is divided into several grid units at certain intervals; temporally, data is collected at fixed time intervals. Multidimensional raw monitoring data refers to data obtained from monitoring the agricultural production area from multiple aspects, covering the physiological state of crops and various factors of the surrounding environment. The crop physiological response time-series data layer records the changes in the physiological state of crops at different times, such as the changes in leaf water potential and stem diameter over time; the environmental factor dynamic data layer reflects the dynamic changes of environmental factors over time, such as changes in soil moisture content and relative humidity.

[0014] In practical implementation, for the time series sequence of crop leaf water potential in the crop physiological response time series data layer, the pressure chamber method can be used to collect data. The pressure chamber method involves placing crop leaves in a sealed pressure chamber and gradually increasing the pressure until water begins to seep from the cut surface; the recorded pressure value is the leaf water potential. For the time series sequence of stem diameter changes, a laser displacement sensor can be used. This sensor emits a laser beam onto the stem surface and measures minute changes in stem diameter based on the time and angle of the reflected light. For the spatial distribution sequence of soil moisture content in the environmental factor dynamic data layer, a time domain reflectometer (TDR) can be used. The TDR calculates soil moisture content based on the propagation speed of the electromagnetic waves in the soil by emitting high-frequency electromagnetic waves into the soil. For the sequence of relative humidity changes, a humidity sensor can be used for real-time monitoring. Humidity sensors are typically based on capacitance or resistance principles, measuring humidity by the change in capacitance or resistance caused by changes in air humidity. By organizing and arranging the collected data according to the time and spatial information of a preset spatiotemporal sampling grid, a monitoring data field with spatiotemporal continuity can be generated.

[0015] Step S200: Extract features from the monitoring data field, extract physiological rhythm features and environmental fluctuation features under different time dimensions, and obtain a crop water-sensitive feature set through feature correlation analysis. The crop water-sensitive feature set contains key feature indicators that characterize changes in crop water status.

[0016] Feature extraction is the process of mining key information reflecting crop physiological characteristics and environmental changes from monitoring data fields. Physiological rhythm characteristics refer to the regularity of physiological activities exhibited by crops at different time scales, such as diurnal variations and medium- to long-term trend changes; environmental fluctuation characteristics refer to the fluctuations of environmental factors at different frequency ranges, including random fluctuations and trend fluctuations. Feature correlation analysis studies the degree of correlation between physiological rhythm characteristics and environmental fluctuation characteristics, identifying feature combinations that have a significant impact on changes in crop water status. The crop water-sensitive feature set is a collection of key feature indicators reflecting changes in crop water status.

[0017] Optionally, step S200 may include the following steps S210 to S250:

[0018] Step S210: Perform spatiotemporal layering analysis on the monitoring data field. Based on the timestamp of data acquisition and spatial sampling coordinate information, decompose the monitoring data field into an independent crop physiological response time series data layer and an environmental factor dynamic data layer. The crop physiological response time series data layer includes the crop leaf water potential time series and the stem diameter change time series. The environmental factor dynamic data layer includes the soil moisture spatial distribution sequence and the air relative humidity change sequence.

[0019] Spatiotemporal hierarchical analysis separates the monitoring data field based on the temporal and spatial information of the data, forming independent data layers for crop physiological response data and environmental factor data. Timestamps record the information at each data acquisition moment, while spatial sampling coordinates indicate the specific location of data collection. The crop leaf water potential time series is a sequence showing the change of crop leaf water potential over time, reflecting the moisture status of crop leaves; the stem diameter change time series records the changes in stem diameter at different times, which is also closely related to the crop's moisture status. The soil moisture content spatial distribution sequence shows the distribution of soil moisture content in different locations within the agricultural production area, while the relative humidity change sequence reflects the dynamic changes in air humidity over time.

[0020] In practice, the timestamp and spatial sampling coordinates of each data point are first extracted from the monitoring data field. Then, the data are classified based on this information: data related to crop physiological state are grouped into the crop physiological response time-series data layer, and data related to environmental factors are grouped into the environmental factor dynamic data layer.

[0021] Step S220: Perform time series decomposition processing on the crop physiological response time series data layer. Use the ensemble empirical mode decomposition method to decompose the crop leaf water potential time series into multiple intrinsic mode function components. Extract the periodic features and amplitude features of each component as physiological rhythm features. The physiological rhythm features include diurnal variation periodic features and medium- and long-term trend features.

[0022] Time series decomposition processing breaks down time series data in crop physiological response time series data layers into different components to better analyze their inherent patterns. Ensemble empirical mode decomposition (EMD) is an adaptive signal decomposition method that can decompose complex time series into a series of intrinsic mode function (IMF) components with different frequencies and amplitudes. IMF components are signal components with specific physical meaning, each representing a particular fluctuation pattern in the time series. Periodicity characteristics refer to the fluctuation period of each component, while amplitude characteristics represent the fluctuation range of each component. The diurnal variation periodicity reflects the daily variation pattern of crop water status, while the medium- to long-term trend characteristics reflect the changing trend of crop water status over a longer period.

[0023] Optionally, step S220 may include the following steps S221 to S224:

[0024] Step S221: Perform extreme point identification processing on the time series of crop leaf water potential, mark the local maxima and local minima in the series, and use cubic spline interpolation to fit the upper and lower envelopes. Calculate the mean of the upper and lower envelopes as the average envelope.

[0025] Extreme point identification involves finding local maximum and minimum values ​​in a crop leaf water potential time series. A local maximum is the maximum value within a local region of the series, and a local minimum is the minimum value within that same region. Cubic spline interpolation is a commonly used method that fits data points using a cubic polynomial in each small interval, resulting in a smoother interpolation curve. The upper envelope connects all local maxima, the lower envelope connects all local minima, and the average envelope is the average of the upper and lower envelopes.

[0026] In practice, the time series of crop leaf water potential is first traversed, comparing the magnitude of each data point with its neighboring data points to identify and mark local maxima and minima. Then, cubic spline interpolation is used to fit the upper envelope using the local maxima as control points and the lower envelope using the local minima as control points. Specifically, the cubic spline interpolation method constructs a cubic polynomial for each inter-cell interval, ensuring that the polynomial's value equals the data point's value at the endpoints of the inter-cell intervals and satisfies certain smoothness conditions. Finally, the average value of the upper and lower envelopes at each time point is calculated to obtain the average envelope.

[0027] Step S222: Subtract the average envelope from the crop leaf water potential time series to obtain the preliminary decomposition components. Perform an orthogonality test on the preliminary decomposition components. When the orthogonality index of the preliminary decomposition components meets the preset orthogonality threshold, it is determined as the first intrinsic mode function component. Otherwise, repeat the extreme point identification and envelope fitting operations on the preliminary decomposition components.

[0028] The preliminary decomposition components are obtained by subtracting the mean envelope from the crop leaf water potential time series, representing the remaining fluctuations after removing the average trend. The orthogonality test checks whether the preliminary decomposition components meet certain orthogonality conditions; the orthogonality index is a numerical value measuring the degree of orthogonality of the preliminary decomposition components. The preset orthogonality threshold is a pre-defined standard value used to determine whether the preliminary decomposition components meet the requirements. If the orthogonality index of the preliminary decomposition component meets the preset orthogonality threshold, the component is considered an intrinsic mode function component; otherwise, extreme point identification and envelope fitting operations need to be performed again on the preliminary decomposition component until an intrinsic mode function component that meets the conditions is obtained.

[0029] In practice, the difference between the crop leaf water potential time series and the mean envelope is first calculated to obtain preliminary decomposition components. Then, methods such as correlation coefficients are used to calculate the orthogonality index of the preliminary decomposition components. For example, the correlation coefficient between the preliminary decomposition components and other determined intrinsic mode function components can be calculated. If the absolute value of the correlation coefficient is less than a preset orthogonality threshold, the orthogonality requirement is considered met. If the requirement is not met, the extreme point identification and envelope fitting of the preliminary decomposition components are performed again to obtain new preliminary decomposition components, and orthogonality testing is performed again. This process is repeated until the first intrinsic mode function component that meets the conditions is found.

[0030] Step S223: After removing the determined intrinsic mode function components, repeat the above decomposition process on the remaining sequence until the remaining sequence becomes a monotonic sequence or a constant sequence, and obtain a decomposition result set containing multiple intrinsic mode function components and a residual component.

[0031] Removing the identified intrinsic modal function components involves subtracting them from the crop leaf water potential time series to obtain the remaining sequence. Then, the remaining sequence undergoes further operations such as extreme point identification, envelope fitting, and orthogonality testing to obtain new intrinsic modal function components. This process is repeated until the remaining sequence becomes a monotonic sequence (i.e., the data in the sequence are monotonically increasing or decreasing) or a constant sequence (the data values ​​in the sequence remain unchanged). The residual components are the final remaining monotonic or constant sequences, representing the long-term trend in the sequence. Specifically, after identifying the first intrinsic modal function component, it is subtracted from the crop leaf water potential time series to obtain the remaining sequence. Then, steps S221-S222 are repeated on the remaining sequence to obtain the second intrinsic modal function component. This process is continuously repeated, removing each new intrinsic modal function component from the remaining sequence.

[0032] Step S224: Perform power spectrum analysis on each intrinsic mode function component in the decomposition result set, calculate the dominant frequency and corresponding amplitude characteristics of each component, determine the component with the dominant frequency within the preset diurnal frequency range as the diurnal variation cycle characteristic, and determine the residual component with the dominant frequency below the preset low frequency threshold as the medium- and long-term trend characteristic. The amplitude characteristics of the diurnal variation cycle characteristic represent the daily variation amplitude of crop water status, and the slope of the medium- and long-term trend characteristic represents the cumulative change trend of crop water status.

[0033] Power spectrum analysis is a method for analyzing the frequency components of a signal. It transforms the signal from the time domain to the frequency domain, thereby obtaining the power distribution of each frequency component. The dominant frequency is the frequency with the highest power in the signal, and the corresponding amplitude characteristic is the magnitude of the signal's amplitude at that frequency. The preset diurnal frequency range is a pre-defined frequency range corresponding to the diurnal cycle, usually around the daily cycle; the preset low-frequency threshold is a frequency boundary used to distinguish between medium- and long-term trends and short-term fluctuations.

[0034] In practice, Fast Fourier Transform (FFT) is used to perform power spectrum analysis on each intrinsic mode function component in the decomposition result set. The power spectrum of each component is calculated using FFT, and the frequency corresponding to the maximum value in the power spectrum is identified as the dominant frequency of that component. Simultaneously, the amplitude value at that frequency is recorded as the corresponding amplitude characteristic. For components whose dominant frequencies fall within a preset diurnal frequency range, they are identified as diurnal variation periodic characteristics, and their amplitude characteristics reflect the range of changes in crop water status within a day. For residual components whose dominant frequencies are below a preset low-frequency threshold, they are identified as medium- to long-term trend characteristics. By calculating the slope of the residual components, the cumulative trend of crop water status changes can be obtained. For example, in a crop with a growth cycle of several months, by decomposing its leaf water potential time series and performing power spectrum analysis, the diurnal variation periodic characteristics and medium- to long-term trend characteristics can be accurately identified, providing an important basis for understanding the changing patterns of crop water status.

[0035] Step S230: Perform multi-scale fluctuation extraction processing on the dynamic data layer of environmental factors. Decompose the spatial distribution sequence of soil moisture content into fluctuation components of different frequency bands through continuous wavelet transform. Calculate the energy proportion and fluctuation amplitude of the fluctuation components of each frequency band as environmental fluctuation features. Environmental fluctuation features include random fluctuation features and trend fluctuation features.

[0036] Multi-scale fluctuation extraction is the process of extracting fluctuation information at different scales from the dynamic data layer of environmental factors. Continuous wavelet transform is a time-frequency analysis method that decomposes a signal at different time and frequency scales, thereby obtaining the fluctuation components of the signal in different frequency bands. Energy proportion refers to the proportion of energy of each frequency band fluctuation component in the total energy, while fluctuation amplitude is the magnitude of the fluctuation component's amplitude. Random fluctuation characteristics are the irregular, random fluctuation components in environmental factors, while trend fluctuation characteristics are the fluctuation components with a certain trend.

[0037] In practice, the first step is to perform a continuous wavelet transform on the spatial distribution sequence of soil moisture content. The continuous wavelet transform uses a wavelet basis function. By performing translation and scaling operations on the wavelet basis function and convolving it with the spatial distribution sequence of soil moisture content, wavelet coefficients at different scales and locations are obtained. These wavelet coefficients reflect the fluctuations of the signal at different frequency bands and time points. Then, the energy of the fluctuation component in each frequency band is calculated based on the wavelet coefficients. The energy can be calculated by the sum of squares of the wavelet coefficients. The energy percentage of each frequency band fluctuation component is calculated, i.e., the energy of that frequency band divided by the total energy. Simultaneously, the maximum and minimum amplitudes of each frequency band fluctuation component are identified, and their difference is calculated as the fluctuation amplitude. For example, for the spatial distribution sequence of soil moisture content in a large farmland, the continuous wavelet transform can decompose it into fluctuation components in multiple frequency bands, such as high-frequency random fluctuation components and low-frequency trend fluctuation components.

[0038] Step S240: Input physiological rhythm features and environmental fluctuation features into the feature association analysis unit to calculate nonlinear correlation degree. Use the mutual information entropy algorithm to calculate the correlation strength value between each type of physiological rhythm feature and environmental fluctuation feature, and generate a feature correlation degree matrix containing the correlation strength values ​​of feature pairs.

[0039] The feature association analysis unit is a module used to analyze the association between physiological rhythm features and environmental fluctuation features. Nonlinear association calculation takes into account the potentially complex nonlinear relationship between physiological rhythm features and environmental fluctuation features, and uses a specific algorithm to calculate the association strength between them. The mutual information entropy algorithm is an algorithm used to measure the dependency between two random variables. It obtains the mutual information entropy value by calculating the joint entropy and marginal entropy of the two variables; the larger the mutual information entropy value, the stronger the association between the two variables. The feature association matrix is ​​a matrix whose rows and columns correspond to the physiological rhythm feature category and the environmental fluctuation feature category, respectively, and the matrix elements are the association strength values ​​of the corresponding feature pairs.

[0040] Optionally, step S240 may include the following steps S241 to S245:

[0041] Step S241: Convert the diurnal variation cycle feature and the medium-to-long-term trend feature in the physiological rhythm features into continuous numerical sequences, and perform the same numerical sequence conversion process on the random fluctuation feature and the trend fluctuation feature in the environmental fluctuation features to obtain a standardized set of physiological feature sequences and a set of environmental feature sequences.

[0042] Continuous numerical sequences are sequences of consecutive numerical values, facilitating subsequent calculations and analysis. Standardization involves unifying different types of characteristic data to give them the same dimensions and range, enabling better comparison and calculation. Physiological characteristic sequence sets are collections of numerical sequences derived from physiological rhythm characteristics, while environmental characteristic sequence sets are collections of numerical sequences derived from environmental fluctuation characteristics.

[0043] In practice, for the diurnal variation cycle features and medium- to long-term trend features in the physiological rhythm characteristics, they are arranged in chronological order, and the feature value corresponding to each time point is used as an element in the sequence to form a continuous numerical sequence. The same processing is performed for the random fluctuation features and trend fluctuation features in the environmental fluctuation characteristics. Then, the resulting sequences are standardized. Standardization can be achieved using the Z-score standardization method, which involves subtracting the mean of the sequence from each element and then dividing by the standard deviation of the sequence. This processing is performed on all physiological rhythm features and environmental fluctuation features to obtain standardized sets of physiological feature sequences and environmental feature sequences.

[0044] Step S242: Select a physiological feature sequence from the set of physiological feature sequences as the target physiological sequence, and select an environmental feature sequence from the set of environmental feature sequences as the target environmental sequence. Construct a feature pair sequence pair, which includes the target physiological sequence and the target environmental sequence.

[0045] The target physiological sequence is selected from the set of physiological feature sequences for association analysis with environmental features, while the target environmental sequence is selected from the set of environmental feature sequences. A feature-sequence pair is a pair of sequences formed by combining the target physiological sequence and the target environmental sequence; it is used to subsequently calculate the strength of the association between them.

[0046] In practice, the physiological feature sequence set and the environmental feature sequence set are iterated through in a loop. Each physiological feature sequence is selected as the target physiological sequence, and each environmental feature sequence is selected as the target environmental sequence, thus constructing feature pair sequences. Each feature pair sequence represents a combination of physiological and environmental features, providing a specific analytical object for subsequent calculation of the correlation strength between them.

[0047] Step S243: Perform joint probability density estimation on the feature pair sequence pairs, use the kernel density estimation method to calculate the two-dimensional joint probability density function of the target physiological sequence and the target environmental sequence, and calculate the marginal probability density function of the target physiological sequence and the marginal probability density function of the target environmental sequence.

[0048] Joint probability density estimation uses known data to estimate the joint probability density function of two random variables, describing the probability distribution of the two variables taking the same value simultaneously. Kernel density estimation is a nonparametric method that places a kernel function around each data point and then sums all the kernel functions to obtain an estimate of the probability density function. The two-dimensional joint probability density function describes the probability distribution of the target physiological sequence and the target environmental sequence taking the same value simultaneously, while the marginal probability density functions describe the probability distribution of the target physiological sequence and the target environmental sequence taking their respective values.

[0049] In practice, for the target physiological sequence and target environmental sequence in the feature pair sequence, a kernel density estimation method is used to estimate the joint probability density. First, a suitable kernel function, such as the Gaussian kernel function, is selected. For each data point, its contribution to the surrounding area is calculated using the kernel function, with that data point as the center. The contributions of all data points are summed to obtain the estimated value of the two-dimensional joint probability density function. Simultaneously, by integrating the two-dimensional joint probability density function over another variable, the marginal probability density functions of the target physiological sequence and the target environmental sequence can be obtained.

[0050] Step S244: Based on the two-dimensional joint probability density function, the marginal probability density function of the target physiological sequence, and the marginal probability density function of the target environmental sequence, calculate the mutual information entropy value of the feature pair sequence pair using the mutual information entropy calculation method, and use the mutual information entropy value as the association strength value of the feature pair.

[0051] The mutual information entropy calculation method uses a two-dimensional joint probability density function, the marginal probability density function of the target physiological sequence, and the marginal probability density function of the target environmental sequence to calculate the mutual information entropy value between two variables. The mutual information entropy value reflects the degree of dependence between the two variables; the larger the value, the stronger the correlation between the two variables.

[0052] In practice, this integral value is calculated using numerical integration to obtain the mutual information entropy value of the feature pair sequence pair. For example, using the trapezoidal integral method, the integration interval is divided into several smaller intervals. An approximate calculation is performed within each smaller interval using the trapezoidal formula, and then the calculation results across all smaller intervals are summed to obtain the mutual information entropy value. This mutual information entropy value is used as the association strength value of the feature pair to measure the degree of association between the target physiological sequence and the target environmental sequence.

[0053] Step S245: Traverse all combinations of physiological feature sequences and environmental feature sequences to generate a feature correlation matrix containing the correlation strength values ​​of all feature pairs. The row dimension of the feature correlation matrix corresponds to the physiological rhythm feature category, the column dimension corresponds to the environmental fluctuation feature category, and the matrix element value is the correlation strength value of the corresponding feature pair.

[0054] In practice, following steps S242-S244, all combinations in the physiological feature sequence set and the environmental feature sequence set are calculated to obtain the association strength value for each feature pair. Then, these association strength values ​​are arranged according to the physiological rhythm feature category and the environmental fluctuation feature category to form a feature association degree matrix. This feature association degree matrix provides a clear visual representation of the degree of association between different physiological rhythm features and environmental fluctuation features, offering an important basis for subsequent screening of crop water-sensitive features.

[0055] Step S250: Sort the correlation strength values ​​in the feature correlation matrix in descending order, select feature combinations with correlation strength values ​​greater than the preset correlation threshold as the initial sensitive feature candidate set, evaluate the feature importance of the initial sensitive feature candidate set through the recursive feature elimination algorithm, calculate the influence weight value of each candidate feature on crop water state changes, and generate the crop water sensitive feature set by combining the feature indicators with the highest influence weight values ​​according to the preset proportion.

[0056] Descending sorting arranges the correlation strength values ​​in the feature correlation matrix from largest to smallest. The preset correlation threshold is a pre-defined standard value used to filter feature combinations strongly correlated with crop water state changes. The initial sensitive feature candidate set consists of feature combinations with correlation strength values ​​greater than the preset correlation threshold. The recursive feature elimination algorithm is a feature selection algorithm that evaluates the importance of each feature by progressively eliminating unimportant features. The influence weight value measures the degree of influence of each candidate feature on crop water state changes. The preset ratio is the proportion used to determine the size of the final crop water sensitive feature set.

[0057] In practice, the correlation strength values ​​in the feature correlation matrix are first sorted in descending order. A sorting algorithm, such as quicksort, can be used for its high efficiency. Then, feature combinations with correlation strength values ​​greater than a preset correlation threshold are selected to form an initial sensitive feature candidate set. Next, a recursive feature elimination algorithm is used to evaluate the feature importance of the initial sensitive feature candidate set. This algorithm typically begins by training a machine learning model, such as a linear regression model, on all candidate features to obtain the coefficients for each feature. Then, the feature with the smallest absolute coefficient is eliminated, the model is retrained, and this process is repeated until a preset stopping condition is met. After each feature elimination, the importance score of each feature is recorded as an influence weight. Finally, the candidate features are ranked according to their influence weights, and the top-ranked feature combinations with a preset proportion are selected to generate the crop water-sensitive feature set.

[0058] Step S300: Input the crop water-sensitive feature set into the crop water metabolism model, infer the coupling relationship between crop root water absorption rate and soil water conductivity, and generate a water supply and demand balance coefficient that characterizes the dynamic balance state of crop water supply and demand.

[0059] A crop water metabolism model is a model used to describe the processes of water absorption, transport, and consumption in crops. It can simulate the water metabolism of crops based on an input set of crop water-sensitive features. Crop root water absorption rate refers to the speed at which crop roots absorb water from the soil, while soil water conductivity is an indicator reflecting the soil's ability to conduct water. Coupling relationship refers to the mutual influence and interaction between crop root water absorption rate and soil water conductivity. The water supply and demand balance coefficient is an indicator used to measure the balance between crop water supply and demand, with a value ranging from complete imbalance to complete balance.

[0060] Optionally, step S300 may include the following steps S310 to S350:

[0061] Step S310: Align the key feature indicators in the crop water-sensitive feature set according to their feature dimensions, and combine them into a feature vector matrix according to the preset feature input order. The row dimension of the feature vector matrix is ​​the number of feature indicators, and the column dimension is the number of time sampling points.

[0062] Feature dimension alignment ensures that key feature indicators in the crop water-sensitive feature set have the same dimension and range, enabling effective combination and computation. The preset feature input order is a pre-determined arrangement of features, and the feature vector matrix is ​​a matrix composed of crop water-sensitive feature sets, where rows represent different feature indicators and columns represent different time sampling points.

[0063] In practice, the first step is to examine the dimensions and ranges of each key feature indicator in the crop water-sensitive feature set. For feature indicators with different dimensions, interpolation or sampling is performed to ensure they have the same number of time sampling points. Then, according to the preset feature input order, the key feature indicators are arranged sequentially to form a feature vector matrix.

[0064] Step S320: Input the feature vector matrix into the feature preprocessing layer of the crop water metabolism model, and perform nonlinear mapping processing on the feature vector matrix through a multilayer perceptron network to reduce the dimensionality of the feature vector matrix to hidden feature vectors. The hidden feature vectors contain comprehensive representation information of crop water status.

[0065] The feature preprocessing layer is the part of the crop water metabolism model used to preprocess the input data, performing operations such as data cleaning, transformation, and feature extraction. A multilayer perceptron (MLP) network is an artificial neural network consisting of an input layer, hidden layers, and an output layer. It can non-linearly map the input data using a non-linear activation function, thereby extracting complex features from the data. The latent feature vector is a low-dimensional vector obtained after processing by the MLP network, containing comprehensive information about the crop's water status.

[0066] In practice, the feature vector matrix is ​​input into the feature preprocessing layer of the crop water metabolism model. The input layer of the multilayer perceptron network receives the feature vector matrix, and the hidden layers contain multiple neurons. Each neuron processes the input through weighted summation and a nonlinear activation function (such as the ReLU function). The ReLU function is expressed as f(x) = max(0,x), which introduces nonlinear factors and enhances the network's expressive power. Through the processing of multiple hidden layers, the feature vector matrix is ​​gradually reduced in dimensionality, and finally, the hidden feature vector is obtained in the output layer.

[0067] Step S330: Call the coupling relationship inference layer of the crop water metabolism model, associate the hidden feature vector with the real-time soil moisture conductivity data, and dynamically model the temporal change process of crop root water absorption rate through the gated recurrent unit network to generate the root water absorption rate time series curve and the soil moisture conductivity change curve.

[0068] The coupling relationship inference layer is the part of the crop water metabolism model used to infer the coupling relationship between crop root water absorption rate and soil water conductivity. Real-time soil water conductivity data is the soil's ability to conduct water measured at the current moment. The root water absorption rate time-series curve reflects the change of crop root water absorption rate over time, while the soil water conductivity change curve reflects the change of soil water conductivity over time.

[0069] Optionally, step S330 may include the following steps S331 to S335:

[0070] Step S331: Align the latent feature vector with the real-time soil moisture conductivity data on the time axis to ensure that the two types of data are consistent at the time sampling points, and generate an aligned joint input sequence. Each time sampling point of the joint input sequence contains the latent feature vector element and the soil moisture conductivity value at the corresponding time.

[0071] Time axis alignment involves matching the latent feature vectors and real-time soil moisture conductivity data according to time, ensuring they correspond at the same sampling points. The joint input sequence is a sequence composed of the aligned latent feature vectors and real-time soil moisture conductivity data, containing both crop water status information and soil moisture conductivity information.

[0072] In practice, the time sampling points of the hidden feature vector and the real-time soil moisture conductivity data are first checked. If the time sampling points of the two are inconsistent, an interpolation method is used for processing. Then, the hidden feature vector and the real-time soil moisture conductivity data are arranged in chronological order, so that the hidden feature vector elements and soil moisture conductivity values ​​at each time sampling point are combined to form a joint input sequence.

[0073] Step S332: Input the joint input sequence into the input layer of the gated recurrent unit network, and perform feature selection and memory update processing on the joint input sequence through the reset gate and the update gate. The reset gate is used to control the retention ratio of historical state information, and the update gate is used to adjust the influence weight of new input information on the current state.

[0074] The input layer of a gated recurrent unit (GRU) network receives a joint input sequence. The reset gate and the update gate are two important gating mechanisms in GRU. The reset gate determines the degree to which historical state information is retained at the current moment, while the update gate controls the degree to which new input information updates the current state.

[0075] In practice, when the combined input sequence enters the input layer of the gated recurrent unit network, the reset gate and update gate begin to operate. The calculation formula for the reset gate is r. t =σ(W r [h t-1 ,x t ]+b r ), where r t It is the output of the reset gate at time t, σ is the sigmoid function, and W r It is the weight matrix of the reset gate, h t-1 It is the hidden state from the previous moment, x t It is the input at the current moment, b r This is the bias term for the reset gate. The output r of the reset gate... t Used to control the retention ratio of historical state information, when r t When r is close to 1, more historical state information is retained; when r is close to 1, more historical state information is retained. t When the value is close to 0, almost no historical state information is retained. The formula for calculating the update gate is z. t =σ(W z [h t-1 ,x t ]+b z ), where z t It is the output of the update gate at time t, W z It is the weight matrix of the updated gate, b z This is the bias term of the update gate. The output z of the update gate... t Used to adjust the weight of the influence of new input information on the current state, when zt When z approaches 1, the new input information has a relatively small impact on the current state; when z approaches 1... t When the value is close to 0, new input information has a greater impact on updating the current state. By working together with the reset gate and the update gate, feature selection and memory update processing are performed on the joint input sequence, enabling the gated recurrent unit network to better process sequence data.

[0076] Step S333: The gated feature information is stored through the cell state unit of the gated recurrent unit network, and the cell state value at the current moment is calculated. The cell state value integrates historical state information and current input feature information.

[0077] The cell state unit is the part of the gated recurrent unit network used for storing and updating information. After gating through reset and update gates, feature information is passed to the cell state unit. The cell state value is the output of the cell state unit at the current moment, combining historical state information and current input feature information.

[0078] In practice, the formula for calculating cell state values ​​is h. t =(1-z t )h ’ t +z t h t-1 , where h t It is the cell state value at the current moment, z. t It updates the output of the gate, h. ’ t =tanh(W[h t-1 ·r t ,x t ]+b) represents the candidate cell state value, · indicates element-wise multiplication, W is the weight matrix, and b is the bias term. First, based on the output r of the reset gate... t Historical state information h t-1 Filter and then compare with the current input x t The values ​​are input together into a linear transformation and a nonlinear activation function (such as the tanh function) to obtain the candidate cell state value h. ’ t Finally, based on the output z of the update gate... t The candidate cell state value h ’ t and historical status information h t-1 By performing a weighted combination, we obtain the cell state value h at the current moment. t .

[0079] Step S334: Input the cell state value into the output layer of the gated recurrent unit network, selectively output the cell state value through the output gate, generate the root water absorption rate prediction value at the current moment, and arrange the root water absorption rate prediction values ​​of all time sampling points in chronological order to generate the root water absorption rate time series curve.

[0080] The output gate is the part of the gated recurrent unit network (RON) used to control information output. It selectively outputs cell state values, only outputting information related to root water uptake rate. The root water uptake rate prediction is the RNU's prediction of the crop root water uptake rate at the current moment. The root water uptake rate time-series curve is formed by arranging the predicted root water uptake rates from all time sampling points in chronological order, visually demonstrating how the crop root water uptake rate changes over time.

[0081] In practical implementation, the calculation formula for the output gate is o t =σ(W o h t +b o ), where o t W is the output of the output gate at time t. o It is the weight matrix of the output gate, b o This is the bias term of the output gate. The output of the output gate is o. t Used for cell state value h t Perform selective output and combine it with a linear transformation W out o t +b out By combining these values, the predicted root water absorption rate at the current moment can be obtained. Arranging the predicted root water absorption rates from all time sampling points in chronological order generates a time-series curve of root water absorption rate.

[0082] Step S335: Perform time series smoothing on the soil moisture conductivity values ​​in the joint input sequence, use the moving average filtering method to eliminate high-frequency noise interference, and generate a smoothed soil moisture conductivity change curve. The time sampling interval of the soil moisture conductivity change curve is consistent with the root water absorption rate time series curve.

[0083] Time series smoothing is a process that manipulates time series data to reduce noise and fluctuations, resulting in smoother data. Specifically, for the soil moisture conductivity values ​​in the joint input sequence, a moving average filtering method is used. This eliminates high-frequency noise interference in the soil moisture conductivity values, generating a smoothed soil moisture conductivity variation curve. Simultaneously, it ensures that the time sampling interval of the soil moisture conductivity variation curve is consistent with the root water absorption rate time series curve.

[0084] Step S340: Perform a synchronization analysis on the time series curve of root water absorption rate and the change curve of soil moisture conductivity, calculate the cross correlation coefficient and phase difference parameter of the two curves, and construct a dynamic coupling relationship model between root water absorption rate and soil moisture conductivity. The dynamic coupling relationship model takes soil moisture conductivity as the input variable and root water absorption rate as the output variable.

[0085] Synchronicity analysis is the process of studying the temporal consistency and correlation between time-series curves of root water absorption rate and soil moisture conductivity. The cross-correlation coefficient is an indicator that measures the correlation between two curves, reflecting the degree of similarity between them at different time delays; the phase difference parameter represents the relative displacement of the two curves over time. The dynamic coupling model describes the dynamic relationship between root water absorption rate and soil moisture conductivity, and can predict changes in root water absorption rate based on changes in soil moisture conductivity.

[0086] In practice, the first step is to calculate the cross-correlation coefficient between the time-series curve of root water absorption rate and the curve of soil moisture conductivity change. The formula for calculating the cross-correlation coefficient can use the commonly used Pearson correlation coefficient method, which will not be elaborated here. By calculating the cross-correlation coefficient at different time delays k, the k value corresponding to the maximum value is found, which is the phase difference parameter. Then, based on the cross-correlation coefficient and the phase difference parameter, a dynamic coupling relationship model between root water absorption rate and soil moisture conductivity is constructed. A linear regression model can be used, such as y = ax + b, where y is the root water absorption rate, x is the soil moisture conductivity, and a and b are coefficients to be determined. The values ​​of coefficients a and b are estimated using the least squares method to minimize the error between the model's predicted values ​​and the actual values.

[0087] Step S350: Calculate the predicted root water absorption rate under different soil moisture conductivity conditions based on the dynamic coupling relationship model. Compare the predicted root water absorption rate with the water consumption rate of crop transpiration to generate a water supply and demand balance coefficient that characterizes the dynamic balance of crop water supply and demand. The numerical range of the water supply and demand balance coefficient is the interval from complete imbalance to complete balance. The closer the value is to the balance endpoint, the closer the water supply and demand is to the balance state.

[0088] Dynamic coupling models can predict root water absorption rates based on input soil moisture conductivity values. Crop transpiration water consumption rate is the rate at which crops lose water through transpiration, and it is related to the crop's growth status and environmental conditions. The water supply and demand balance coefficient is an indicator that measures the degree of balance between crop water supply (root water absorption rate) and demand (crop transpiration water consumption rate).

[0089] Optionally, step S350 may include the following steps S351 to S354:

[0090] Step S351: Select multiple discrete test points evenly within the preset range of soil moisture conductivity values, input the soil moisture conductivity values ​​of each discrete test point into the dynamic coupling relationship model, obtain the predicted value of root water absorption rate under the corresponding conditions, and generate a table of correspondence between soil moisture conductivity and root water absorption rate.

[0091] The preset range of soil moisture conductivity values ​​is a pre-defined interval of possible values ​​for soil moisture conductivity based on actual conditions. Discrete test points are specific soil moisture conductivity values ​​uniformly selected within this range. The soil moisture conductivity-root water absorption rate correspondence table records the predicted root water absorption rate values ​​corresponding to different soil moisture conductivity values.

[0092] In practice, firstly, a preset range for soil moisture conductivity is determined, for example, from 0.1 to 1.0. Then, multiple discrete test points are evenly selected within this range, such as 10 test points with values ​​of 0.1, 0.2, 0.3, ..., 1.0. The soil moisture conductivity value of each discrete test point is input into a dynamic coupling relationship model, and the predicted root water absorption rate under the corresponding conditions is calculated according to the model's calculation formula. Finally, the soil moisture conductivity values ​​and the corresponding predicted root water absorption rates are compiled into a table to generate a table showing the correspondence between soil moisture conductivity and root water absorption rate.

[0093] Step S352: Collect the crop transpiration rate under the current environmental conditions as the rate of water consumption by crop transpiration per unit land area.

[0094] Crop canopy transpiration rate refers to the rate at which a crop population per unit area loses water through transpiration under current environmental conditions. It is influenced by various factors, including crop species, growth stage, ambient temperature, light intensity, and air humidity.

[0095] In practice, a lysimeter is used to collect the crop canopy transpiration rate. A lysimeter is an instrument that measures soil moisture evaporation and crop transpiration, calculating the crop canopy transpiration rate by measuring the decrease in soil moisture over a certain period. In operation, the lysimeter is installed in the farmland, ensuring its measurement range covers a certain area of ​​the crop canopy. Over a period of time, the change in soil moisture recorded in the lysimeter is recorded, along with the measurement time and area. The crop canopy transpiration rate is calculated using the formula E = ΔW / (A×t), where E is the crop canopy transpiration rate, ΔW is the decrease in soil moisture, A is the measurement area, and t is the measurement time.

[0096] Step S353: Based on the soil moisture conductivity-root water absorption rate correspondence table and real-time soil moisture conductivity measurement, the predicted root water absorption rate under the current soil moisture conditions is calculated using a linear interpolation method. The predicted root water absorption rate represents the rate at which crop roots absorb water from the soil.

[0097] Real-time soil moisture conductivity measurement is data on the soil's ability to conduct water, measured at the current moment. In practice, the interval containing the real-time soil moisture conductivity measurement is first located from the soil moisture conductivity-root water absorption rate correspondence table. Assuming the real-time soil moisture conductivity measurement value x lies between two adjacent soil moisture conductivity values ​​x1 and x2 in the correspondence table, the corresponding predicted root water absorption rates are y1 and y2, respectively. Then, the predicted root water absorption rate y under the current soil moisture conditions is calculated using the linear interpolation formula y = y1 + (y2 - y1)(x - x1) / (x2 - x1).

[0098] Step S354: Compare the predicted root water absorption rate with the crop transpiration water consumption rate to generate a water supply and demand balance coefficient that characterizes the dynamic balance of crop water supply and demand. Determine the value of the water supply and demand balance coefficient based on the comparison results.

[0099] In practice, the difference between the predicted root water absorption rate and the crop transpiration rate is first calculated as Δ = RT, where R is the predicted root water absorption rate and T is the crop transpiration rate. Then, the water supply and demand balance coefficient is determined based on the magnitude of the difference Δ. A balance range can be set; for example, when Δ is between -0.1 and 0.1, water supply and demand are considered to be close to equilibrium, and the water supply and demand balance coefficient can be set between 0.8 and 1.0. When Δ is less than -0.1, it indicates insufficient crop water supply, and the water supply and demand balance coefficient can be set between 0 and 0.2. When Δ is greater than 0.1, it indicates excessive crop water supply, and the water supply and demand balance coefficient can be set between 0.2 and 0.8. The specific coefficient setting can be adjusted according to the actual situation so that the water supply and demand balance coefficient can accurately reflect the dynamic balance state of crop water supply and demand.

[0100] Step S400: Construct an irrigation decision optimization model based on the water supply and demand balance coefficient and real-time soil moisture data. The irrigation decision optimization model includes an irrigation total optimization model that optimizes the total irrigation amount with the goal of maximizing water saving rate and a growth stage irrigation quota allocation model that allocates irrigation quotas for each growth stage with the goal of maximizing water use efficiency. The irrigation parameter set is obtained through the irrigation total optimization model and the growth stage irrigation quota allocation model.

[0101] The water supply and demand balance coefficient reflects the balance between crop water supply and demand, while real-time soil moisture data indicates the current soil moisture content. Irrigation decision optimization models are used to formulate optimal irrigation plans, comprehensively considering both water conservation and water use efficiency. The total irrigation optimization model aims to maximize water conservation by adjusting the total irrigation water volume for the total irrigated area; the growth stage irrigation quota allocation model aims to maximize water use efficiency by rationally allocating irrigation water for each crop growth stage. The irrigation parameter set is a collection of parameters including total irrigation volume and irrigation quotas for each growth stage.

[0102] Optionally, step S400 may include the following steps S410 to S440:

[0103] Step S410: Perform spatiotemporal matching processing on the water supply and demand balance coefficient and the real-time soil moisture data to ensure that the two types of data are consistent in spatial sampling location and temporal sampling point, and generate a matched joint decision input dataset. The joint decision input dataset contains multiple sampling unit data in the spatial dimension.

[0104] Spatiotemporal matching aligns the water supply-demand balance coefficient and real-time soil moisture data in time and space, ensuring they correspond at the same location and time. Spatial sampling location refers to the specific location where data is collected, while temporal sampling point records the time of data collection. The joint decision input dataset consists of the matched water supply-demand balance coefficient and real-time soil moisture data, containing data from multiple sampling units, each corresponding to a specific spatial location.

[0105] In practice, the spatial and temporal sampling locations of the water supply-demand balance coefficient and real-time soil moisture data are first checked. If the spatial sampling locations are inconsistent, spatial interpolation methods, such as Kriging interpolation, are used to estimate data at unknown locations based on data from known sampling points. If the temporal sampling points are inconsistent, temporal interpolation methods, such as linear interpolation, are used. Then, the matched water supply-demand balance coefficient and real-time soil moisture data are arranged according to their spatial and temporal sampling locations to form a joint decision input dataset.

[0106] Step S420: Construct an irrigation total quantity optimization model with the goal of maximizing water saving rate. The decision variable of the irrigation total quantity optimization model is the total irrigation water volume of the total irrigated area. The constraints include the upper limit constraint of soil moisture, the lower limit constraint of crop water demand, and the flow constraint of the irrigation system. Solve the irrigation total quantity optimization model through a genetic algorithm to obtain the optimal irrigation total quantity value.

[0107] Water-saving rate refers to the ratio of water saved through optimized irrigation schemes to traditional irrigation water volume. Maximizing the water-saving rate means reducing irrigation water volume as much as possible while meeting crop growth needs. The decision variable is the total irrigation water volume required to be determined in the total irrigation volume optimization model; here, it is the total irrigation water volume for the total irrigated area. The upper limit constraint on soil moisture stipulates that the soil moisture content after irrigation cannot exceed the field capacity to prevent problems such as root hypoxia caused by excessive soil moisture. The lower limit constraint on crop water demand requires that the soil moisture content in the crop root zone after irrigation not be lower than the wilting coefficient to ensure normal crop growth. The irrigation system flow constraint limits the irrigation water volume per unit time to not exceed the maximum water supply capacity of the irrigation system.

[0108] Optionally, step S420 may include the following steps S421 to S426:

[0109] Step S421: Define the objective function of the total irrigation volume optimization model as maximizing the water-saving rate. The water-saving rate is calculated based on the comparison between the traditional total irrigation volume and the optimized total irrigation volume. The traditional total irrigation volume is the conventional irrigation water volume based on experience, and the optimized total irrigation volume is the model decision variable.

[0110] The objective function is the function value that needs to be maximized or minimized in the total irrigation volume optimization model; in this case, it is the water-saving rate. The formula for calculating the water-saving rate is S = (V c -V o ) / V c Where S is the water saving rate, and V c It is the total amount of traditional irrigation, V o It involves optimizing the total irrigation volume. Traditional irrigation volume is determined based on past experience and habits, while optimizing the total irrigation volume is achieved by solving a model to obtain the optimal irrigation volume.

[0111] Step S422: Set constraints for the total irrigation optimization model. The upper limit constraint on soil moisture stipulates that the soil moisture content after irrigation shall not exceed the field capacity. The lower limit constraint on crop water requirement stipulates that the soil moisture content in the crop root zone after irrigation shall not be lower than the wilting coefficient. The flow constraint of the irrigation system stipulates that the irrigation water volume per unit time shall not exceed the maximum water supply capacity of the irrigation system.

[0112] Field holding capacity refers to the maximum water content that soil can retain after drainage under gravity; it represents the upper limit of soil moisture. The wilting index refers to the soil moisture content at which crops begin to permanently wilt due to water shortage; it represents the lower limit of crop water requirements. The maximum water supply capacity of an irrigation system is the maximum amount of water that the system can provide per unit time.

[0113] In practical implementation, the upper limit constraint on soil moisture is determined based on the soil type and characteristics, specifically the field capacity. For example, for a certain loam soil, its field capacity is 30% (volume water content). When optimizing the total irrigation volume, it is ensured that the soil moisture content after irrigation does not exceed 30%. The lower limit constraint on crop water requirements is determined based on the crop type and growth stage. For example, if the wilting coefficient of a crop in its mid-growth stage is 10%, then the soil moisture content in the crop root zone after irrigation must not be lower than 10%. For irrigation system flow constraints, the specifications and performance of the irrigation system are understood to determine its maximum water supply capacity. By setting these constraints, it is ensured that the irrigation scheme meets the crop's growth needs without placing an excessive burden on the soil and irrigation system.

[0114] Step S423: Initialize the parameter settings of the genetic algorithm, including population size, number of generations, crossover probability and mutation probability. Randomly generate the initial population. Each individual corresponds to a candidate solution for total irrigation amount. Individual encoding adopts real number encoding method.

[0115] Population size refers to the number of individuals in each generation of the genetic algorithm. The number of generations is the number of iterations the algorithm performs. Crossover probability is the probability of a crossover operation occurring during genetic operations, and mutation probability is the probability of a mutation operation occurring. The initial population is a set of candidate solutions for total irrigation amount randomly generated at the start of the algorithm; each candidate solution represents a possible irrigation amount scheme. Real-number encoding means that the total irrigation amount for each individual is directly represented by a real number.

[0116] In practice, the population size and number of generations are determined based on the complexity of the problem and computational resources. When randomly generating the initial population, a random number generator is used to generate the value for each individual, based on the range of the total irrigation amount.

[0117] Step S424: Evaluate the fitness of each individual in the initial population, use the objective function value as the fitness function value, standardize the fitness value, and use the roulette wheel selection method to select individuals with high fitness as parent individuals.

[0118] Fitness assessment is the process of evaluating the merits of each individual in a population. The fitness function value reflects the degree to which an individual is adapted to a target. In the total irrigation optimization model, the value of the objective function (water-saving rate) is used as the fitness function value; that is, the higher the water-saving rate, the higher the fitness of the individual. Standardization involves uniformly processing fitness values ​​to give them the same dimensions and range, facilitating comparison and selection. The roulette wheel selection method is a probability-based selection method where the probability of each individual being selected is proportional to its fitness value.

[0119] In practice, for each individual in the initial population, its corresponding total irrigation amount is substituted into the objective function (water-saving rate calculation formula) to calculate the fitness function value. Then, the fitness values ​​of all individuals are standardized. One standardization method is to divide each individual's fitness value by the sum of all individual fitness values ​​to obtain each individual's relative fitness value. When using the roulette wheel selection method, each individual's relative fitness value is considered as the area of ​​a sector on the roulette wheel; the larger the fitness value, the larger the sector area. By randomly spinning the roulette wheel, the individual that lands on the sector is selected as the parent individual.

[0120] Step S425: Use the single-point crossover operator to perform crossover operation on the selected parent individuals to generate offspring individuals, and use the Gaussian mutation operator to perform mutation operation on the offspring individuals to maintain population diversity.

[0121] The single-point crossover operator is a genetic operation that randomly selects a crossover point on the chromosome of a parent individual and exchanges the portions of the chromosomes of two parents after the crossover point, thus generating two offspring individuals. The Gaussian mutation operator is a mutation operation that randomly selects one or more gene loci on the chromosome of an offspring individual and mutates the gene values ​​at these loci, with the mutation magnitude following a Gaussian distribution. Maintaining population diversity is to avoid the algorithm getting trapped in local optima and to ensure that the algorithm can search a wider solution space.

[0122] In practice, for the selected parent individuals, a single-point crossover operator is used for crossover. First, a crossover point is randomly selected. For the generated offspring individuals, a Gaussian mutation operator is used for mutation. One or more gene loci are randomly selected on the chromosome of the offspring individual, and a random mutation amplitude is generated according to a Gaussian distribution to adjust the gene value at that locus. Through this mutation operation, the gene values ​​of the offspring individuals change, increasing the diversity of the population. By continuously performing crossover and mutation operations, new offspring individuals are generated, driving the evolution of the population.

[0123] Step S426: Repeat the selection, crossover and mutation operations until the preset number of generations is reached, and decode the individual with the highest fitness value in the last generation population as the optimal total irrigation value.

[0124] Repeated selection, crossover, and mutation operations are the core iterative process of genetic algorithms. In each generation, superior individuals are selected based on their fitness values ​​for reproduction, and new individuals are generated through crossover and mutation, continuously optimizing the quality of the population. A preset number of generations is the termination condition for the algorithm's iteration; when this number of generations is reached, the algorithm stops iterating. The final generation is the last generation of the population, where the individual with the highest fitness value represents the currently found optimal solution.

[0125] In practice, starting from the initial population, selection, crossover, and mutation operations are performed according to steps S424-S425 to generate a new generation of population. The fitness of each individual in the new generation is evaluated, and then selection, crossover, and mutation operations are performed again. This process is repeated until the preset number of generations is reached.

[0126] In the last generation of the population, identify the individual with the highest fitness value. Since individuals are encoded using real numbers, this individual's value is directly used as the optimal total irrigation value.

[0127] Step S430: Construct an irrigation quota allocation model for each growth stage with the optimization objective of maximizing water use efficiency. The decision variables of the irrigation quota allocation model for each growth stage are the irrigation water allocation ratios for each growth stage of the crop. The constraints include water sensitivity index constraints and total irrigation constraints for each growth stage. Solve the irrigation quota allocation model for each growth stage using the particle swarm optimization algorithm to obtain the irrigation quota ratios for each growth stage.

[0128] Water use efficiency refers to the yield generated by a crop per unit of water consumed. Maximizing water use efficiency means increasing crop yield as much as possible given a certain amount of irrigation water. The decision variable is the irrigation quota allocation ratio for each growth stage that needs to be determined in the growth stage model, i.e., the proportion of irrigation water allocated to each growth stage of the crop. The water sensitivity index constraint for each growth stage stipulates that the proportion of irrigation water allocated to different growth stages must be positively correlated with the water sensitivity index of that stage. The total irrigation constraint requires that the sum of the irrigation water for each growth stage equals the optimal total irrigation value.

[0129] Optionally, step S430 may include the following steps S431-S436:

[0130] Step S431: Define the objective function of the irrigation quota allocation model for each growth stage as maximizing water use efficiency. The water use efficiency is calculated as the ratio of crop yield to total irrigation water volume. Crop yield is calculated from the irrigation water volume of each growth stage through the crop water production function model.

[0131] The objective function is the value that needs to be maximized in the irrigation quota allocation model for different growth stages, i.e., water use efficiency. The crop water production function model describes the relationship between crop yield and irrigation water volume at each growth stage, and can predict crop yield based on irrigation water volume at different growth stages.

[0132] In practical implementation, the first step is to determine the crop water production function model. This model can be derived from extensive experimental data and statistical analysis, taking into account the crop's water requirements and responses at different growth stages. For example, by conducting experiments on the yield of a certain crop under different irrigation volumes, a function model containing multiple parameters can be established. This model can calculate the crop yield based on the irrigation volume at each growth stage.

[0133] The crop yield is calculated by substituting the irrigation water allocation ratios for each growth stage into the crop water production function model. Then, according to the method for calculating water use efficiency, the crop yield is divided by the total irrigation water volume to obtain the water use efficiency value. This water use efficiency value is used as the objective function value, with the goal of maximizing water use efficiency by adjusting the irrigation water allocation ratios for each growth stage.

[0134] Step S432: Set constraints for the irrigation quota allocation model for each growth stage. The water sensitivity index constraint for each growth stage stipulates that the irrigation water allocation ratio for different growth stages must be positively correlated with the water sensitivity index of that stage. The total irrigation amount constraint stipulates that the sum of irrigation water for each growth stage is equal to the optimal total irrigation amount.

[0135] The water sensitivity index at each growth stage reflects the crop's sensitivity to water at different growth stages. The higher the water sensitivity index, the more sensitive the crop is to water demand at that stage, requiring a larger allocation of irrigation water. The total irrigation quota constraint is to ensure that the total irrigation water volume does not exceed the optimal total irrigation value obtained through the total irrigation quota optimization model when allocating irrigation quotas for each growth stage.

[0136] In practice, the water sensitivity index for each growth stage is determined based on the crop type and growth characteristics. For example, a crop may have a lower water sensitivity index during the seedling stage and a higher water sensitivity index during the flowering stage. When allocating irrigation water, the proportion of irrigation water allocated during the flowering stage should be relatively high, and the proportion allocated during the seedling stage should be relatively low, to meet the water requirements of the crop at different growth stages.

[0137] For total irrigation constraints, after determining the irrigation water allocation ratio for each growth stage, each ratio is multiplied by the optimal total irrigation value to obtain the actual irrigation water volume for each growth stage. Then, it is checked whether the sum of the actual irrigation water volumes for each growth stage equals the optimal total irrigation value. If they are not equal, the irrigation water allocation ratio for each growth stage needs to be adjusted until the total irrigation constraint condition is met.

[0138] Step S433: Initialize the parameter settings of the particle swarm optimization algorithm, including particle swarm size, learning factor, inertia weight and maximum number of iterations. Each particle corresponds to an irrigation quota allocation ratio vector, and the vector elements are the allocation ratios for each growth stage.

[0139] The particle swarm size refers to the number of particles in the swarm. The learning factor is a parameter that controls the particle's flight towards its historical best position and the global best position. The inertia weight is a parameter that controls the particle's flight speed. The maximum number of iterations is the termination condition for the algorithm's iteration. Each particle represents a possible irrigation quota allocation scheme, and the irrigation quota allocation ratio vector is a vector composed of the irrigation water allocation ratio for each growth stage.

[0140] In practice, the particle swarm size and maximum number of iterations are determined based on the complexity of the problem and computational resources. An irrigation quota allocation vector is randomly initialized for each particle. Each vector element takes a value between 0 and 1, and the sum of all elements equals 1.

[0141] Step S434: Perform a feasibility check on the position vector of each particle to ensure that all constraints are met. Repair infeasible solutions by using the penalty function method to penalize particles that violate the constraints and reduce their fitness values.

[0142] Feasibility checks examine whether the irrigation quota allocation scheme represented by each particle satisfies the constraints of the irrigation quota allocation model for each growth stage, including water sensitivity index constraints and total irrigation volume constraints. Infeasible solutions are those that do not meet the constraints and require remediation to make them satisfy the conditions. The penalty function method is a way to punish particles that violate the constraints by reducing their fitness value, thus decreasing the probability of these particles being selected in subsequent iterations.

[0143] In practice, for each particle's position vector, the first step is to check if the total irrigation amount constraint is met. Each element in the vector is multiplied by the optimal total irrigation amount value to calculate the sum of the actual irrigation water volume for each growth stage, and then compared to the optimal total irrigation amount value. If they are not equal, the vector elements are adjusted to meet the total irrigation amount constraint. For example, the vector elements can be scaled proportionally so that the sum of the elements equals 1, and multiplying by the optimal total irrigation amount value satisfies the total irrigation amount requirement. Next, the constraint of water sensitivity index for each growth stage is checked. Based on the water sensitivity index for each growth stage, the rationality of the irrigation water allocation ratio is determined. If it is not rational, the vector elements are adjusted so that the irrigation water allocation ratio is positively correlated with the water sensitivity index. For example, the irrigation water allocation ratio for growth stages with high water sensitivity indices is increased, and the irrigation water allocation ratio for growth stages with low water sensitivity indices is decreased.

[0144] For particles that still violate the constraints after inspection, a penalty function method is used to impose a penalty. The penalty function can be set according to the degree of constraint violation; the greater the degree of violation, the heavier the penalty. For example, for a particle that violates the total irrigation amount constraint, the difference between its actual irrigation amount and the optimal total irrigation amount is multiplied by a penalty coefficient, and this penalty value is subtracted from its fitness value, thus reducing its fitness value. In this way, the particle is encouraged to fly towards a feasible solution that satisfies the constraints.

[0145] Step S435: Calculate the fitness value of each particle, use the objective function value as the fitness function value, track the individual optimal position and global optimal position of the particle, and iteratively update the velocity and position of the particle according to the velocity update formula and the position update formula.

[0146] The fitness value reflects the quality of the irrigation quota allocation scheme represented by each particle. The objective function value (water use efficiency) is used as the fitness function value; a higher fitness value indicates a better scheme. The individual optimal position is the best position reached by the particle itself historically, while the global optimal position is the best position reached by all particles in the entire particle swarm historically. The velocity update formula and position update formula are used in particle swarm optimization algorithms to update the particle's flight velocity and position.

[0147] In practice, for each particle, its position vector (irrigation quota allocation ratio vector) is substituted into the objective function of the irrigation quota allocation model for the reproductive stage to calculate its fitness value. Then, the particle's current fitness value is compared with the fitness value of its historical best position. If the current fitness value is higher, the individual's best position is updated to the current position. At the same time, the fitness values ​​of all particles in the entire particle swarm are compared, and the particle with the highest fitness value is identified, with its position being taken as the global best position.

[0148] The velocity and position of each particle are updated using the velocity update formula and the position update formula. The velocity update formula considers the particle's current velocity, its tendency to fly towards its historical best position, and its tendency to fly towards the global best position. The position update formula updates the particle's position based on the updated velocity. For example, in the current iteration, a particle calculates its new velocity using the velocity update formula, and then adds the new velocity to its current position using the position update formula to obtain its new position. By continuously iterating and updating the particle's velocity and position, the particle swarm gradually converges towards the global best position.

[0149] Step S436: When the maximum number of iterations is reached or the number of consecutive preset times of the global optimal position remains unchanged, stop the iteration and use the allocation ratio vector corresponding to the global optimal position as the irrigation quota ratio for each growth stage.

[0150] The maximum number of iterations is a pre-defined termination condition for the algorithm. When the algorithm reaches this number of iterations, it stops iterating. Another termination condition is that the global optimal position remains unchanged for a pre-defined number of consecutive iterations. This indicates that the global optimal position has not changed in a number of consecutive iterations, meaning the algorithm has converged to a stable solution.

[0151] In practice, after each iteration, it is checked whether the maximum number of iterations has been reached. If the maximum number of iterations has been reached, iteration stops. Simultaneously, the changes in the global optimal position are recorded, and the number of times the global optimal position remains unchanged is counted. When this number reaches a preset number, iteration also stops. After the algorithm stops iterating, the allocation ratio vector corresponding to the global optimal position is used as the irrigation quota ratio for each growth stage. Through iterative optimization using the particle swarm optimization algorithm, an irrigation quota allocation scheme that maximizes water use efficiency for each growth stage under the given constraints is obtained.

[0152] Step S440: Multiply the optimal total irrigation value with the irrigation quota ratio of each growth stage to obtain the specific irrigation quota value of each growth stage. Combine the spatial zoning information to generate an irrigation parameter set containing the total irrigation amount, the growth stage quota, and the zoning irrigation ratio. The element values ​​of the irrigation parameter set are non-negative real numbers.

[0153] After obtaining the optimal total irrigation value and the irrigation quota ratio for each growth stage, multiplying the two yields the specific irrigation quota value for each growth stage. Spatial zoning information refers to the division of agricultural production areas into different regions based on factors such as soil type, topography, and crop planting conditions. Each region may have different irrigation needs. Zoning irrigation ratio refers to the proportion of irrigation water in each zone to the total irrigation water.

[0154] In practice, the optimal total irrigation amount is first multiplied by the irrigation quota ratio for each growth stage to obtain the specific irrigation quota value for each growth stage. Combining spatial zoning information, the zoning irrigation ratio is determined based on factors such as soil conditions and crop growth in each zone. For example, if the agricultural production area is divided into three zones, the zoning irrigation ratios for zones one, two, and three are 0.3, 0.4, and 0.3, respectively. The specific irrigation quota values ​​for each growth stage are then allocated according to the zoning irrigation ratios to obtain the specific irrigation water volume for each zone at each growth stage. Finally, the total irrigation amount, the specific irrigation quota values ​​for each growth stage, and the zoning irrigation ratios are combined to generate an irrigation parameter set.

[0155] Step S500: Drive the irrigation actuator to perform zonal variable irrigation operations based on the irrigation parameter set, collect crop transpiration rate change data after irrigation, correct the water supply and demand balance coefficient based on the crop transpiration rate change data, and generate the final irrigation control command to drive the irrigation actuator to perform zonal variable irrigation operations.

[0156] The irrigation parameter set includes information such as total irrigation volume, growth stage quotas, and zonal irrigation ratios. This information drives irrigation execution agencies to perform zonal variable irrigation operations. Zonal variable irrigation refers to adjusting irrigation volume and timing based on the irrigation needs of different zones. Crop transpiration rate changes reflect changes in crop water status after irrigation. Analyzing this data can correct the water supply and demand balance coefficient, making it more accurately reflect the actual water requirements of crops. The final irrigation control instructions are generated based on the corrected water supply and demand balance coefficient and are used to guide irrigation execution agencies in subsequent irrigation operations.

[0157] Optionally, step S500 may include the following steps S510-S560:

[0158] Step S510: Analyze the zonal irrigation ratio information in the irrigation parameter set, divide the agricultural production area into multiple irrigation control zones, and each zone corresponds to an independent irrigation execution mechanism control unit.

[0159] Zonal irrigation ratio information refers to the information on the allocation of irrigation water to different zones within the irrigation parameter set. By analyzing this information, agricultural production areas can be divided into multiple irrigation control zones according to certain rules. Each irrigation control zone has similar soil conditions, crop growth, and irrigation needs, and each zone corresponds to an independent irrigation execution mechanism control unit, facilitating independent control of irrigation in each zone.

[0160] In practice, the first step is to extract the zoning irrigation ratio information from the irrigation parameter set. Then, based on factors such as the topography, soil type, and crop distribution of the agricultural production area, and combined with the zoning irrigation ratio information, the agricultural production area is divided into multiple irrigation control zones. For example, a farmland may be divided into three irrigation control zones based on differences in soil moisture: Zone 1 has low soil moisture, Zone 2 has moderate soil moisture, and Zone 3 has high soil moisture. Each zone is equipped with an independent irrigation actuator control unit, such as a solenoid valve or water pump controller, through which the irrigation water volume and irrigation time of each zone can be independently controlled.

[0161] Step S520: Calculate the irrigation amount per unit area based on the irrigation quota value and area of ​​each zone, and convert the irrigation amount per unit area into control parameters of the irrigation actuator. The control parameters include irrigation time length and irrigation flow rate.

[0162] The irrigation quota for each zone is the amount of irrigation water required for that zone, determined based on the set of irrigation parameters and the zone's irrigation ratio. The zone area is the actual area of ​​each irrigation control zone. Irrigation volume per unit area refers to the amount of water required to irrigate each unit area of ​​land, obtained by dividing the irrigation quota for each zone by its area. The control parameters for the irrigation actuators are used to control their operation, including irrigation duration and flow rate.

[0163] Step S530: Send control commands to the irrigation actuators of each zone in sequence according to the preset irrigation order, drive the irrigation actuators to perform the zone variable irrigation operation, and record the actual irrigation start time and end time of each zone.

[0164] The preset irrigation sequence is a pre-determined order in which different zones will be irrigated, which can be set according to factors such as crop growth needs and soil moisture conditions. Control commands are instructions containing information such as irrigation duration and flow rate, used to control the operation of the irrigation actuators. Recording the actual start and end times of irrigation for each zone facilitates subsequent evaluation and analysis of the irrigation operation's effectiveness.

[0165] In practice, according to the preset irrigation sequence, starting with the first zone, control commands are sent to the irrigation actuators in that zone. These commands are transmitted to the control unit of the irrigation actuator via a communication interface, such as through a wireless communication module to the solenoid valve or water pump controller. The irrigation actuator then begins operation based on the control commands. Simultaneously, the actual irrigation start time for that zone is recorded.

[0166] Once the irrigation operation in one zone is completed, the actual irrigation end time is recorded. Then, following the preset irrigation sequence, control commands are sent to the irrigation actuators of the second zone, and the above process is repeated until the irrigation operation in all zones is completed. In this way, the orderly execution of irrigation operations in different zones is achieved.

[0167] Step S540: Within a preset time interval after the irrigation operation is completed, collect crop transpiration rate data for each zone. Set multiple sampling points for each zone. The collection time covers a complete photoperiod to generate a crop transpiration rate change dataset. The crop transpiration rate change dataset contains transpiration rate measurements over time.

[0168] The preset time interval refers to a period of waiting after irrigation to allow the crop sufficient time to respond to the irrigation before collecting crop transpiration rate data. Setting multiple sampling points in each zone improves data representativeness and ensures accurate reflection of crop transpiration within that zone. A complete photoperiod refers to the entire period from sunrise to sunset; collecting data covering a complete photoperiod provides a comprehensive understanding of transpiration rate changes under different light conditions. The crop transpiration rate variation dataset consists of crop transpiration rate measurements from each zone at different time points, reflecting changes in crop water status after irrigation.

[0169] Optionally, step S540 may include the following steps S541-S545:

[0170] Step S541: Set up multiple sampling points in each irrigation control zone according to the principle of uniform distribution. Select representative crop plants with consistent growth status as the measurement object for each sampling point, and record the geographical coordinates and crop growth stage information of the sampling point.

[0171] The principle of uniform distribution means that sampling points are evenly distributed throughout each irrigation control zone to ensure comprehensive coverage of the crop conditions within that zone. Representative crop plants are those with good growth and exhibiting common characteristics of crops in that zone; selecting such plants as measurement subjects improves the accuracy of the results. Recording the geographical coordinates of the sampling points and the crop growth stage information provides a reference for subsequent data analysis and comparison.

[0172] In practice, for each irrigation control zone, the number and location of sampling points are determined based on the zone's area and shape. For example, for a large zone, 10 sampling points can be set up and evenly distributed within the zone. Representative crop plants with similar growth characteristics are selected at each sampling point, such as similar plant height, number of leaves, and color. The geographical coordinates of the sampling points are recorded using a Global Positioning System (GPS) device, while simultaneously observing the crop's growth stages, such as seedling, flowering, and fruiting stages, and recording relevant information.

[0173] Step S542: After the irrigation operation is completed, the crop transpiration rate data will be collected at a preset time period. A portable photosynthesis instrument will be used to continuously measure at each sampling point. Data will be collected at preset time intervals, from the preset time before sunrise to the preset time after sunset, covering the complete photoperiod.

[0174] The preset time period is a waiting period after irrigation to ensure the crop responds significantly to the irrigation during this time. A portable photosynthesis meter is an instrument used to measure crop photosynthetic and transpiration rates, enabling rapid and accurate measurement of crop physiological parameters. The preset time interval refers to the time interval between each data collection session. By continuously measuring according to this interval, the change in crop transpiration rate over time can be obtained. Starting from the preset duration before sunrise and ending at the preset duration after sunset, covering the complete photoperiod, a comprehensive understanding of the changes in crop transpiration rate under different light conditions can be achieved.

[0175] Step S543: At each measurement time point, the transpiration rate of multiple crop plants at each sampling point is measured, and the arithmetic mean of multiple measurements is taken as the representative value of the transpiration rate of the sampling point at that time point. The ambient temperature and light intensity data at the time of measurement are recorded.

[0176] Measuring the transpiration rate of multiple crop plants at each sampling point can reduce measurement errors and improve the accuracy of the results. Taking the arithmetic mean of multiple measurements as the representative value of the transpiration rate at that sampling point and time point can more accurately reflect the crop transpiration situation at that sampling point. Recording the ambient temperature and light intensity data during the measurement allows for analysis of the impact of environmental factors on crop transpiration rate.

[0177] In practice, transpiration rates were measured on multiple crop plants at each sampling point at each measurement time point. Ambient temperature and light intensity data were recorded using a thermometer and a light intensity sensor. These data, along with representative transpiration rate values, were recorded to provide more information for subsequent data analysis.

[0178] Step S544: Arrange the representative values ​​of evaporation rate of each sampling point at different time points in chronological order to generate a time series sequence of evaporation rate for a single sampling point. Average the time series sequences of evaporation rate of all sampling points in the same partition to obtain the partition average evaporation rate time series sequence.

[0179] Arranging the representative transpiration rate values ​​of each sampling point at different time points in chronological order clearly shows the change of crop transpiration rate over time at each sampling point. Averaging the transpiration rate time series of all sampling points within the same region yields the average transpiration rate time series of that region, which more accurately reflects the crop transpiration situation in that region.

[0180] In practice, for each sampling point, the representative values ​​of its evaporation rate at different time points are arranged in chronological order to form a time series.

[0181] Step S545: Combine the time series of average transpiration rates of all partitions into a crop transpiration rate change dataset. The row dimension of the crop transpiration rate change dataset is the number of partitions, the column dimension is the number of time sampling points, and the data element value is the average transpiration rate measurement of the corresponding partition at the corresponding time point.

[0182] By combining the time-series average transpiration rates of all partitions, a complete crop transpiration rate variation dataset can be formed, facilitating comprehensive analysis of crop transpiration across all partitions. The row dimension of the crop transpiration rate variation dataset represents the number of partitions, the column dimension represents the number of time sampling points, and the data element value is the average transpiration rate measurement for each partition at each time point.

[0183] In practice, the time series of the average transpiration rate for each zone is arranged sequentially to form a matrix. The number of rows in the matrix equals the number of zones, and the number of columns equals the number of time sampling points. For example, if there are 3 irrigation control zones, and the time series of the average transpiration rate for each zone is 24 (i.e., 24 time sampling points), then the crop transpiration rate variation dataset is a 3×24 matrix, where each element represents the average transpiration rate measurement for the corresponding zone at the corresponding time point. This dataset allows for a direct observation of the changes in crop transpiration rate across different zones throughout the entire photoperiod.

[0184] Step S550: Perform trend analysis on the crop transpiration rate change dataset, calculate the transpiration rate change rate before and after irrigation, and adjust the water supply and demand balance coefficient accordingly according to the comparison results of the transpiration rate change rate and the preset threshold based on the preset rules.

[0185] Trend analysis involves analyzing datasets of crop transpiration rate changes to understand the trend of transpiration rate changes over time. The rate of change in transpiration rate before and after irrigation reflects the degree of change in crop water status after irrigation. The preset threshold is a pre-defined standard value used to determine whether a change in transpiration rate is significant. The preset rule is a rule used to adjust the water supply and demand balance coefficient based on the comparison between the rate of change in transpiration rate and the preset threshold.

[0186] Optionally, step S550 may include the following steps S551-S554:

[0187] Step S551: Extract the transpiration rate data of the same period before irrigation from the crop transpiration rate change dataset as a reference dataset, and align the transpiration rate data after irrigation with the reference dataset on the time axis to ensure that the time sampling points of the two sets of data correspond one-to-one.

[0188] Pre-irrigation transpiration rate data refers to the crop's transpiration rate before irrigation, under the same time period and light conditions as after irrigation. Aligning the post-irrigation transpiration rate data with a reference dataset along the time axis ensures that the two sets of data are compared at the same point in time, improving the accuracy of the comparison results.

[0189] In practice, the transpiration rate data for the same period before irrigation is extracted from the crop transpiration rate variation dataset and used as a reference dataset. For example, if the data after irrigation is the transpiration rate data from 1 hour before sunrise to 1 hour after sunset, then the transpiration rate data for the same time period before irrigation is extracted from historical data and used as a reference dataset.

[0190] The transpiration rate data after irrigation and the reference dataset are aligned on the time axis. If the time sampling points of the two sets of data are not exactly the same, interpolation methods, such as linear interpolation, are used to ensure that the two sets of data have corresponding values ​​at the same time points. Time axis alignment ensures that the two sets of data can be accurately compared.

[0191] Step S552: Calculate the change in evaporation rate at each time sampling point, that is, subtract the evaporation rate measurement value at the corresponding time point in the reference dataset from the evaporation rate measurement value after irrigation, and generate a sequence of evaporation rate changes.

[0192] The change in transpiration rate reflects the change in the crop's transpiration rate at each time point after irrigation relative to before irrigation. By calculating the change in transpiration rate at each time point, a sequence of transpiration rate changes can be obtained, which shows the trend of crop transpiration rate throughout the entire photoperiod.

[0193] In practice, for each time sampling point, the evaporation rate measurement value after irrigation is subtracted from the corresponding time point evaporation rate measurement value in the reference dataset to obtain the change in evaporation rate at that time point. The changes in evaporation rate at all time points are then arranged in chronological order to generate a sequence of evaporation rate changes.

[0194] Step S553: ​​Integrate the transpiration rate change sequence to calculate the total transpiration rate change within the complete photoperiod. Compare the total transpiration rate change with the total transpiration rate of the reference dataset to obtain the transpiration rate change rate.

[0195] Integration involves summing the transpiration rate change sequence over the entire photoperiod to obtain the total transpiration rate change over the complete photoperiod. The total transpiration rate in the reference dataset is the sum of the transpiration rate measurements at all time points in the reference dataset. The transpiration rate change rate is the ratio of the total transpiration rate change to the total transpiration rate in the reference dataset, reflecting the degree of change in crop transpiration rate after irrigation relative to before irrigation.

[0196] In practice, numerical integration methods are used to integrate the transpiration rate change sequence. For example, the trapezoidal integration method is used to divide the transpiration rate change sequence into multiple small intervals. The trapezoidal formula is used to approximate the calculation in each small interval, and then the calculation results of all intervals are added together to obtain the total transpiration rate change over the complete photoperiod.

[0197] The total transpiration rate of the reference dataset is calculated by summing the transpiration rate measurements at all time points in the reference dataset. Then, the total change in transpiration rate is divided by the total transpiration rate of the reference dataset to obtain the rate of change in transpiration rate.

[0198] Step S554: Compare the rate of change of transpiration rate with a preset threshold. The preset threshold is determined based on the crop type and growth stage. When the rate of change of transpiration rate exceeds the preset threshold, it is determined that the current water supply and demand balance coefficient underestimates the crop's water demand, and the water supply and demand balance coefficient is adjusted upward according to the preset rules. When the rate of change of transpiration rate does not exceed the preset threshold, it is determined that the current water supply and demand balance coefficient overestimates the crop's water demand, and the water supply and demand balance coefficient is adjusted downward according to the preset rules. The adjustment range is determined based on the magnitude of the rate of change of transpiration rate, so that the adjusted water supply and demand balance coefficient can more accurately reflect the actual water supply and demand status of the crop.

[0199] The preset threshold is a standard value determined based on the crop type and growth stage. Different crops and growth stages have different water requirements, therefore the preset threshold also varies. When the rate of change in transpiration rate exceeds the preset threshold, it indicates a significant change in the crop's transpiration rate after irrigation, suggesting that the current water supply and demand balance coefficient may not adequately account for the crop's water needs, and the coefficient needs to be adjusted upwards. Conversely, when the rate of change in transpiration rate does not exceed the preset threshold, it indicates a relatively small change in the crop's transpiration rate after irrigation, suggesting that the current water supply and demand balance coefficient may overestimate the crop's water needs, and the coefficient needs to be adjusted downwards. The adjustment range is determined by the magnitude of the rate of change in transpiration rate; the larger the rate of change, the larger the adjustment.

[0200] Step S560: Regenerate the irrigation parameter set based on the adjusted water supply and demand balance coefficient, convert it into the final irrigation control command, which includes the irrigation time, flow rate and duration parameters of each zone, and send the final irrigation control command to the control unit of the irrigation actuator to drive the actuator to perform the zone variable irrigation operation.

[0201] The adjusted water supply and demand balance coefficient more accurately reflects the actual water requirements of crops. Regenerating the irrigation parameter set based on this coefficient allows for a more rational irrigation plan. The final irrigation control command converts the regenerated irrigation parameter set into specific control instructions, including parameters such as irrigation time, flow rate, and duration for each zone. These parameters directly control the operation of the irrigation actuators.

[0202] In practice, the adjusted water supply and demand balance coefficient is substituted into the irrigation decision optimization model to re-optimize the total irrigation volume and allocate irrigation quotas for each growth stage, resulting in a new set of irrigation parameters. The new set of irrigation parameters includes information such as the total irrigation volume that better reflects the actual water needs of crops, irrigation quotas for each growth stage, and the proportion of irrigation in different zones.

[0203] The new set of irrigation parameters is converted into final irrigation control instructions. Based on the irrigation quota and area of ​​each zone, the irrigation volume per unit area of ​​each zone is calculated. Then, combined with the performance of the irrigation equipment, the irrigation time, flow rate, and duration parameters of each zone are determined.

[0204] The final irrigation control command is sent to the control unit of the irrigation actuator via a communication interface, such as through a wireless communication module to the solenoid valve or water pump controller. The irrigation actuator then performs zoned variable irrigation operations based on the final irrigation control command, achieving precise control of crop water, improving irrigation efficiency, and achieving deep water conservation and control.

[0205] It is understood that the various algorithms involved in the embodiments of this invention can all be obtained from relevant content in the prior art. To save space, they will not be elaborated on in the embodiments of this invention. In addition, those skilled in the art can supplement the details based on common knowledge in the art when implementing the solution of this invention. For example, they can use normalization to eliminate dimensional conflicts before feature fusion, use interpolation to eliminate dimensional differences, reasonably set thresholds based on historical data, experience or business scenario requirements, train the model based on a general model training method, set the number of layers in the model structure based on actual needs, select activation functions, etc. This invention will not provide redundant descriptions of overly detailed implementation processes.

[0206] In one embodiment, a computer system, which may be a server, is provided, and its internal structure is shown in Figure 3. The computer system includes a processor, memory, input / output interfaces (I / O), and a communication interface. The processor, memory, and I / O interfaces are connected via a system bus, and the communication interface is connected to the system bus via the I / O interfaces. The processor provides computing and control capabilities. The memory includes non-volatile storage media and internal memory. The non-volatile storage media stores the operating system, computer programs, and a database. The internal memory provides an environment for the operation of the operating system and computer programs stored in the non-volatile storage media. The I / O interfaces are used for exchanging information between the processor and external devices. The communication interface is used for communication with external monitoring devices via a network connection. When the computer program is executed by the processor, it implements a method for deep water-saving and water-controlling agricultural irrigation based on multi-source data fusion.

[0207] Those skilled in the art will understand that the structure shown in Figure 3 is merely a block diagram of a portion of the structure related to the present invention and does not constitute a limitation on the computer system to which the present invention is applied. A specific computer system may include more or fewer components than those shown in the figure, or combine certain components, or have different component arrangements.

[0208] In one embodiment, a computer system is also provided, including a memory and a processor, wherein the memory stores a computer program, and the processor executes the computer program to implement the steps in the above method embodiments.

Claims

1. A method for deep water-saving and water-controlling in agricultural irrigation based on multi-source data fusion, characterized in that, The method includes: collecting multi-dimensional raw monitoring data from agricultural production areas to generate a monitoring data field with spatiotemporal continuity, wherein the monitoring data field includes a crop physiological response time-series data layer and an environmental factor dynamic data layer; extracting features from the monitoring data field, extracting physiological rhythm features and environmental fluctuation features under different time dimensions, and obtaining a crop water-sensitive feature set through feature correlation analysis, wherein the crop water-sensitive feature set includes key feature indicators characterizing changes in crop water status; inputting the crop water-sensitive feature set into a crop water metabolism model, inferring the coupling relationship between crop root water absorption rate and soil water conductivity, and generating a water supply and demand balance coefficient characterizing the dynamic balance state of crop water supply and demand; based on the... An irrigation decision optimization model is constructed using water supply and demand balance coefficients and real-time soil moisture data. This model includes an irrigation total volume optimization model that optimizes the total irrigation volume with the goal of maximizing water-saving rate, and a growth stage irrigation quota allocation model that allocates irrigation quotas for each growth stage with the goal of maximizing water use efficiency. An irrigation parameter set is obtained through these two models. Based on this set of irrigation parameters, the irrigation execution agency is driven to perform zoned variable irrigation operations. Data on crop transpiration rate changes after irrigation are collected. The water supply and demand balance coefficient is corrected based on this crop transpiration rate change data, and a final irrigation control command is generated to drive the irrigation execution agency to perform zoned variable irrigation operations.

2. The method according to claim 1, characterized in that, The monitoring data field undergoes feature extraction, extracting physiological rhythm features and environmental fluctuation features across different time dimensions. A crop water-sensitive feature set is obtained through feature correlation analysis. This set includes key indicator features characterizing changes in crop water status, including: spatiotemporal hierarchical analysis of the monitoring data field, decomposing it into independent crop physiological response time-series data layers and environmental factor dynamic data layers based on data acquisition timestamps and spatial sampling coordinates. The crop physiological response time-series data layer includes crop leaf water potential time-series sequences and stem diameter change time-series sequences, while the environmental factor dynamic data layer includes soil moisture spatial distribution sequences and relative humidity change sequences. The crop physiological response time-series data layer undergoes time-series decomposition, decomposing the crop leaf water potential time-series sequence into multiple intrinsic modal function components. The periodic and amplitude features of each component are extracted as physiological rhythm features, including diurnal variations. The system extracts periodic and medium-to-long-term trend characteristics. It performs multi-scale fluctuation extraction on the dynamic data layer of environmental factors, decomposing the spatial distribution sequence of soil moisture content into fluctuation components of different frequency bands. The energy proportion and fluctuation amplitude of each frequency band fluctuation component are calculated as environmental fluctuation characteristics, including random fluctuation characteristics and trend fluctuation characteristics. The physiological rhythm characteristics and environmental fluctuation characteristics are input into a feature correlation analysis unit for nonlinear correlation degree calculation. The correlation strength value between each type of physiological rhythm characteristic and environmental fluctuation characteristic is calculated, generating a feature correlation degree matrix containing the correlation strength values ​​of feature pairs. Based on the correlation strength values ​​in the feature correlation degree matrix, the system sorts them in descending order, selecting feature combinations with correlation strength values ​​greater than a preset correlation threshold as the initial sensitive feature candidate set. The initial sensitive feature candidate set is evaluated for feature importance, and the influence weight value of each candidate feature on crop water state changes is calculated. The feature indicators with the highest influence weight values ​​(previously a preset proportion) are combined to generate a crop water sensitive feature set.

3. The method according to claim 2, characterized in that, The process of performing time series decomposition on the crop physiological response time series data layer decomposes the crop leaf water potential time series into multiple intrinsic mode function components, extracting the periodic and amplitude features of each component as physiological rhythm features. This includes: identifying extreme points in the crop leaf water potential time series, marking local maxima and minima in the sequence, fitting upper and lower envelopes, and calculating the mean of the upper and lower envelopes as the average envelope; subtracting the average envelope from the crop leaf water potential time series to obtain preliminary decomposed components; performing orthogonality tests on the preliminary decomposed components; and determining the first intrinsic mode function component as the orthogonality index of the preliminary decomposed component if it meets a preset orthogonality threshold. Otherwise, the remaining components are not considered as intrinsic mode functions. The initial decomposition components are repeatedly subjected to extreme point identification and envelope fitting operations. After removing the identified intrinsic mode function components, the above decomposition process is repeated on the remaining sequence until the remaining sequence becomes a monotonic sequence or a constant sequence, resulting in a decomposition result set containing multiple intrinsic mode function components and a residual component. Power spectrum analysis is performed on each intrinsic mode function component in the decomposition result set to calculate the dominant frequency and corresponding amplitude characteristics of each component. Components with dominant frequencies within a preset diurnal frequency range are identified as diurnal variation periodic characteristics, and residual components with dominant frequencies below a preset low-frequency threshold are identified as medium- to long-term trend characteristics. The amplitude characteristics of the diurnal variation periodic characteristics represent the daily variation amplitude of crop water status, and the slope of the medium- to long-term trend characteristics represents the cumulative variation trend of crop water status.

4. The method according to claim 2, characterized in that, The process of inputting the physiological rhythm features and environmental fluctuation features into a feature association analysis unit for nonlinear correlation degree calculation, calculating the correlation strength value between each type of physiological rhythm feature and environmental fluctuation feature, and generating a feature correlation degree matrix containing the correlation strength values ​​of feature pairs includes: converting the diurnal variation periodic features and medium-to-long-term trend features in the physiological rhythm features into continuous numerical sequences, and performing the same numerical sequence conversion process on the random fluctuation features and trend fluctuation features in the environmental fluctuation features to obtain a standardized set of physiological feature sequences and a set of environmental feature sequences; selecting a physiological feature sequence from the set of physiological feature sequences as a target physiological sequence, and selecting an environmental feature sequence from the set of environmental feature sequences as a target environmental sequence, constructing a feature pair sequence pair, wherein the feature pair sequence pair contains the target physiological sequence. The system analyzes the physiological sequence and the target environment sequence. It performs joint probability density estimation on the feature pair sequences, calculates the two-dimensional joint probability density function of the target physiological sequence and the target environment sequence, and calculates the marginal probability density functions of the target physiological sequence and the target environment sequence. Based on the two-dimensional joint probability density function, the marginal probability density function of the target physiological sequence, and the marginal probability density function of the target environment sequence, it calculates the mutual information entropy value of the feature pair sequence, and uses the mutual information entropy value as the association strength value of the feature pair. It then iterates through all combinations of physiological feature sequences and environmental feature sequences to generate a feature correlation matrix containing the association strength values ​​of all feature pairs. The row dimension of the feature correlation matrix corresponds to the physiological rhythm feature category, the column dimension corresponds to the environmental fluctuation feature category, and the matrix element values ​​are the association strength values ​​of the corresponding feature pairs.

5. The method according to claim 1, characterized in that, The process of inputting the crop water-sensitive feature set into a crop water metabolism model to infer the coupling relationship between crop root water absorption rate and soil moisture conduction rate, and generating a water supply and demand balance coefficient characterizing the dynamic balance state of crop water supply and demand, includes: aligning the feature dimensions of each key feature indicator in the crop water-sensitive feature set, combining them into a feature vector matrix according to a preset feature input order, wherein the row dimension of the feature vector matrix is ​​the number of feature indicators, and the column dimension is the number of time sampling points; inputting the feature vector matrix into the feature preprocessing layer of the crop water metabolism model, performing nonlinear mapping processing on the feature vector matrix through a multilayer perceptron network to reduce the dimensionality of the feature vector matrix to latent feature vectors, wherein the latent feature vectors contain comprehensive characterization information of crop water status; and calling the coupling relationship inference layer of the crop water metabolism model to perform correlation modeling between the latent feature vectors and real-time soil moisture conduction rate data, through a gated loop single... Meta-networks dynamically model the temporal variation of crop root water absorption rate, generating a time-series curve of root water absorption rate and a curve of soil moisture conductivity variation. Synchronization analysis is performed on these two curves, calculating their cross-correlation coefficient and phase difference parameter to construct a dynamic coupling model between root water absorption rate and soil moisture conductivity. This model uses soil moisture conductivity as the input variable and root water absorption rate as the output variable. Based on this model, predicted root water absorption rates are calculated under different soil moisture conductivity conditions. These predicted rates are then compared with the rate of water consumption by crop transpiration, generating a water supply-demand balance coefficient that characterizes the dynamic equilibrium of crop water supply and demand. The coefficient ranges from complete imbalance to complete equilibrium; the closer the value is to the equilibrium endpoint, the closer the water supply and demand are to a balanced state.

6. The method according to claim 5, characterized in that, The process of invoking the coupling relationship inference layer of the crop water metabolism model to associate the hidden feature vector with real-time soil moisture conductivity data and dynamically modeling the temporal change process of crop root water absorption rate through a gated recurrent unit network to generate a root water absorption rate time-series curve and a soil moisture conductivity change curve includes: aligning the hidden feature vector with the real-time soil moisture conductivity data along the time axis to generate an aligned joint input sequence, where each time sampling point of the joint input sequence contains the hidden feature vector element and the soil moisture conductivity value at the corresponding time; and inputting the joint input sequence into the input layer of the gated recurrent unit network, where feature selection and memory update processing are performed on the joint input sequence through reset and update gates. The gated feature information is stored in the cell state unit of the gated recurrent unit network, and the cell state value at the current moment is calculated. The cell state value integrates historical state information and current input feature information. The cell state value is input into the output layer of the gated recurrent unit network, and the cell state value is selectively output through the output gate to generate the root water absorption rate prediction value at the current moment. The root water absorption rate prediction values ​​of all time sampling points are arranged in chronological order to generate the root water absorption rate time series curve. The soil moisture conductivity values ​​in the joint input sequence are subjected to time series smoothing to generate a smoothed soil moisture conductivity change curve. The time sampling interval of the soil moisture conductivity change curve is consistent with that of the root water absorption rate time series curve.

7. The method according to claim 5, characterized in that, The process of calculating the predicted root water absorption rate under different soil moisture conductivity conditions based on the dynamic coupling relationship model, comparing the predicted root water absorption rate with the crop transpiration water consumption rate, and generating a water supply and demand balance coefficient characterizing the dynamic balance of crop water supply and demand includes: uniformly selecting multiple discrete test points within a preset range of soil moisture conductivity values, inputting the soil moisture conductivity values ​​of each discrete test point into the dynamic coupling relationship model to obtain the predicted root water absorption rate under the corresponding conditions, and generating a soil moisture conductivity-root water absorption rate correspondence table; collecting current environmental data... The crop canopy transpiration rate under the given conditions is used as the rate at which crop water is consumed by transpiration per unit land area. Based on the soil moisture conductivity-root water absorption rate correspondence table and real-time soil moisture conductivity measurements, a predicted root water absorption rate is calculated under the current soil moisture conditions. This predicted root water absorption rate represents the rate at which crop roots absorb water from the soil. The predicted root water absorption rate is compared with the rate at which crop water is consumed by transpiration to generate a water supply and demand balance coefficient that characterizes the dynamic balance of crop water supply and demand. The value of the water supply and demand balance coefficient is determined by the comparison results.

8. The method according to claim 1, characterized in that, The process involves constructing an irrigation decision optimization model based on the water supply and demand balance coefficient and real-time soil moisture data. This model, along with the irrigation total quantity optimization model and the irrigation quota allocation model for different growth stages, yields a set of irrigation parameters. The process includes: performing spatiotemporal matching processing on the water supply and demand balance coefficient and real-time soil moisture data to ensure consistency between the two types of data in spatial sampling locations and temporal sampling points, generating a matched joint decision input dataset containing multiple sampling units in the spatial dimension; constructing an irrigation total quantity optimization model with maximizing water saving rate as the optimization objective. The decision variable of the irrigation total quantity optimization model is the total irrigation water volume for the total irrigated area, with constraints including upper limit constraints on soil moisture, lower limit constraints on crop water demand, and irrigation system flow constraints. A genetic algorithm is then used to optimize the data. The optimal total irrigation amount is obtained by solving the total irrigation amount optimization model. A growth stage irrigation quota allocation model is constructed with the goal of maximizing water use efficiency. The decision variables of this model are the irrigation water allocation ratios for each growth stage of the crop, and the constraints include water sensitivity index constraints and total irrigation amount constraints for each growth stage. The model is solved using a particle swarm optimization algorithm to obtain the irrigation quota ratios for each growth stage. The optimal total irrigation amount is multiplied by the irrigation quota ratios for each growth stage to obtain the specific irrigation quota values ​​for each growth stage. An irrigation parameter set containing total irrigation amount, growth stage quotas, and zonal irrigation ratios is generated by combining spatial zoning information. The elements of this irrigation parameter set are non-negative real numbers.

9. The method according to claim 1, characterized in that, The process of driving the irrigation actuator to perform zonal variable irrigation operations based on the irrigation parameter set, collecting crop transpiration rate change data after irrigation, correcting the water supply and demand balance coefficient based on the crop transpiration rate change data, and generating a final irrigation control command to drive the irrigation actuator to perform zonal variable irrigation operations includes: parsing the zonal irrigation ratio information in the irrigation parameter set, dividing the agricultural production area into multiple irrigation control zones, each zone corresponding to an independent irrigation actuator control unit; calculating the irrigation amount per unit area based on the irrigation quota value and zone area of ​​each zone, converting the irrigation amount per unit area into control parameters for the irrigation actuator, the control parameters including irrigation duration and irrigation flow rate; sending control commands to the irrigation actuators of each zone sequentially according to a preset irrigation sequence to drive the irrigation actuators to perform zonal variable irrigation operations, and recording the actual irrigation start time of each zone. The start and end times are determined. Within a preset time interval after the irrigation operation is completed, crop transpiration rate data for each zone are collected. Multiple sampling points are set for each zone, and the collection time covers a complete photoperiod, generating a crop transpiration rate change dataset containing time-series transpiration rate measurements. Trend analysis is performed on the crop transpiration rate change dataset to calculate the rate of change of transpiration rate before and after irrigation. Based on the comparison between the rate of change of transpiration rate and a preset threshold, the water supply and demand balance coefficient is adjusted accordingly according to preset rules. The irrigation parameter set is regenerated based on the adjusted water supply and demand balance coefficient and converted into a final irrigation control command. The final irrigation control command contains irrigation time, flow rate, and duration parameters for each zone. The final irrigation control command is sent to the control unit of the irrigation actuator to drive the actuator to perform zoned variable irrigation operations.

10. A computer system, characterized in that, include: processor; And a memory, wherein the memory stores a computer program that, when run by the processor, causes the processor to perform the method as described in any one of claims 1 to 9.

Images