A big data-based intelligent agricultural control system

Abstract

The application relates to the technical field of intelligent agricultural control, and particularly discloses an intelligent agricultural control system based on big data, which collects soil conductivity time series data, yield data, elevation data, meteorological data and variety cultivation attributes in multiple growth seasons; generates a whole field intrinsic conductivity distribution map after removing high-frequency disturbance of conductivity; separates variety effects and cultivation effects based on causal inference to obtain a corrected yield distribution map; superimposes the intrinsic conductivity, the corrected yield and the elevation, and generates a fine management area boundary through compression embedding and spatial constraint clustering; takes a benchmark reference area as a reference, determines unique fertilizer amount and irrigation threshold values of each management area according to a comprehensive coefficient, and generates a variable operation prescription map; the application eliminates the interference of short-term environmental fluctuations and human management on management area division, reduces the uneven growth phenomenon, and realizes water and fertilizer distribution according to needs.

Description

Technical Field

[0001] This invention relates to the field of smart agriculture control technology, and more specifically to a smart agriculture control system based on big data. Background Technology

[0002] Currently, in the field of smart agriculture, there are technical solutions that utilize geographic information systems to overlay multi-year yield distribution maps, soil electrical conductivity maps, and elevation maps, and generate fine-grained management zones within fields using unsupervised clustering algorithms (such as K-means clustering and fuzzy C-means clustering). These solutions typically assume that soil electrical conductivity can stably characterize the intrinsic soil fertility, and that multi-year yield data only reflects spatial differences in soil fertility.

[0003] The existing technology has the following shortcomings: In the complex environment of paddy fields in the south, existing technologies cannot simultaneously remove short-term environmental fluctuations (caused by irrigation, drainage, and rainfall) in the soil electrical conductivity time series signal, nor can they separate the human management effects brought about by variety and farming methods from multi-year yield data. This leads to unsupervised clustering algorithms misjudging high-frequency disturbances and differences in variety and farming as permanent differences in soil fertility, thus generating unreasonable fine management zones. Ultimately, this results in uneven growth of rice in different zones within the same field after variable fertilization. Summary of the Invention

[0004] The purpose of this invention is to provide a smart agricultural control system based on big data to solve the problems mentioned above.

[0005] The objective of this invention can be achieved through the following technical solutions: A smart agriculture control system based on big data includes: The data acquisition module is used to collect time-series data of soil electrical conductivity for multiple consecutive growing seasons, historical yield data, elevation data, meteorological data for each grid cell, as well as rice variety attributes and cultivation method attributes corresponding to each growing season. The electrical conductivity denoising module separates the soil electrical conductivity time series data into multiple fluctuation components, removes the high-frequency disturbance components and retains the low-frequency trend components, and generates an intrinsic electrical conductivity distribution map of the entire field by spatial interpolation of the low-frequency trend components. The yield correction module uses the historical yield data of each grid cell as the dependent variable, rice variety attributes and cultivation method attributes as processing variables, and intrinsic conductivity distribution map, elevation data and meteorological data as covariates. It separates the variety effect and cultivation effect through causal inference, and subtracts the variety effect and cultivation effect from the historical yield data to obtain the corrected yield distribution map. The management area division module overlays the intrinsic conductivity distribution map, the corrected yield distribution map, and the elevation data to form a multidimensional feature for each grid. It then compresses the multidimensional feature into a low-dimensional embedding space, simultaneously optimizes clustering allocation and feature reconstruction within the low-dimensional embedding space, and applies embedding similarity constraints to adjacent grids in the space to generate a fine-grained management area boundary. The prescription map generation module determines the unique fertilization amount and irrigation threshold for each management zone based on the median corrected yield and median intrinsic conductivity within each intensive management zone, and generates a variable operation prescription map.

[0006] As a further aspect of the present invention: the generation of the intrinsic conductivity distribution map of the entire field specifically includes: Adaptive time-frequency decomposition was performed on the soil electrical conductivity time series data of each sampling point to obtain multiple decomposition components arranged from low to high frequency; Calculate the average period of each decomposition component, identify the decomposition components with an average period of less than two hours as high-frequency disturbance components and remove them, and retain the decomposition components with an average period of more than 24 hours as low-frequency trend components. Using the unbiased interpolation method with optimal spatial autocorrelation, the low-frequency trend components of all sampling points are interpolated to each grid cell of the entire field to generate an intrinsic conductivity distribution map.

[0007] As a further aspect of the present invention: the adaptive time-frequency decomposition of the soil conductivity time-series data at each sampling point specifically includes: Identify all local extrema in the conductivity time series data, connect the local maxima and local minima using cubic spline interpolation to form upper and lower envelopes, and calculate the mean curves of the upper and lower envelopes; The mean curve is subtracted successively from the original conductivity time series data. Each subtraction yields a decomposed component. The envelope and mean calculation process is repeated until the remaining signal is monotonic or has no obvious fluctuations, thus obtaining multiple decomposed components. Count the number of times each decomposed component crosses zero within a unit of time, and arrange the decomposed components in order of the number of zero crossings from least to most, thus obtaining the decomposed components arranged from low to high frequency.

[0008] As a further aspect of the present invention: obtaining the corrected production distribution map specifically includes: Based on the intrinsic conductivity distribution, elevation, and meteorological data of each grid cell, the probability value of each grid cell receiving the current rice variety attribute treatment is calculated, and grid cells with similar probability values ​​are grouped into the same comparison group. Within the same control group, the average difference in yield was calculated for different variety attributes and different tillage attributes, and the average difference was used as the estimate of variety effect and tillage effect. By subtracting the estimated effects corresponding to the actual variety attributes and farming methods of the current year from the historical yield data of each grid cell, a corrected yield distribution map is obtained.

[0009] As a further aspect of the present invention: the calculation of the probability value of each grid cell receiving the current rice variety attribute processing specifically includes: The intrinsic conductivity distribution, elevation, and meteorological values ​​of each grid cell are combined to form the three-dimensional feature points of the grid cell; Calculate the standardized Euclidean distance between the three-dimensional feature point and the three-dimensional feature points of all other grid cells in the field, and select a preset number of neighboring grid cells in ascending order of distance; The proportion of neighboring grid cells that share the same rice variety attribute as the current grid cell is used as the probability value for the current grid cell to receive the current variety attribute. Sort all grid cells according to their probability values ​​from low to high, and group adjacent grid cells whose probability values ​​differ from a preset threshold into the same comparison group.

[0010] As a further aspect of the present invention: the generation of the fine-grained management area boundary specifically includes: A three-layer nonlinear compression transformation is constructed to sequentially map the multidimensional features of each grid into intermediate hidden layer representations and the final low-dimensional embedding points; In each iteration, the soft assignment probability from each low-dimensional embedding point to each candidate cluster center is calculated. At the same time, each low-dimensional embedding point is reconstructed back to the original multi-dimensional features, and the similarity between the features before and after reconstruction is used as a constraint to adjust the parameters of the compression transformation. For spatially adjacent grids, calculate the Euclidean distance between their low-dimensional embedding points, and then multiply the Euclidean distance by the spatial weight coefficient and add it to the clustering loss. When the change in soft assignment probability between two consecutive iterations is less than a preset threshold, the boundary of the fine management area to which each grid belongs is output.

[0011] As a further aspect of the present invention: the parameter for adjusting the compression transformation by using the similarity between features before and after reconstruction as a constraint specifically includes: Each low-dimensional embedding point is sequentially subjected to a three-layer expansion transformation that is mirror-symmetric to the compression transformation, gradually restoring the dimension until a reconstructed feature with the same dimension as the original multidimensional feature is obtained. Calculate the absolute value of the difference between the original multidimensional features and the reconstructed features in each dimension, and then take the reciprocal of the weighted sum of the absolute values ​​of the differences in each dimension to obtain the reconstructed similarity. With the goal of maximizing the reconstruction similarity, the projection direction and bias of each layer of the compression transformation are corrected in reverse, so that the reconstruction similarity increases monotonically in each iteration.

[0012] As a further aspect of the present invention: the generation of variable job prescription diagrams specifically includes: The corrected median yield of all finely managed areas in the field is sorted from high to low, and the management area with the highest value is selected as the benchmark reference area. For each non-reference area, calculate the first ratio of the corrected median output of the fine management area to the corrected median output of the reference area, and the second ratio of the median intrinsic conductivity of the fine management area to the median intrinsic conductivity of the reference area. The arithmetic mean of the first ratio and the second ratio is used as the comprehensive coefficient. The standard fertilizer application rate of the benchmark reference area is multiplied by the comprehensive coefficient to obtain the fertilizer application rate of the intensive management area. The benchmark irrigation threshold of the benchmark reference area is subtracted from the product of the comprehensive coefficient and the benchmark irrigation threshold to obtain the irrigation threshold of the intensive management area. The location boundaries of each fine management zone are associated with the corresponding fertilization amount and irrigation threshold to generate a variable operation prescription map.

[0013] The beneficial effects of this invention are: (1) By removing high-frequency disturbance components in the soil electrical conductivity time series data and correcting the yield differences caused by variety and cultivation method, the generated fine management area can truly reflect the soil's intrinsic fertility. This avoids the traditional method from misjudging the differences in variety or cultivation as differences in soil fertility due to short-term environmental fluctuations and human management interference, thereby reducing the uneven growth phenomenon that occurs after farmers apply fertilizer based on incorrect prescription maps.

[0014] (2) Using the median proportion of the benchmark reference area as a comprehensive coefficient, the amount of fertilizer and the irrigation threshold of each management area are determined in a differentiated manner, so that areas with low soil fertility automatically receive less fertilizer and a lower irrigation start threshold, while areas with high soil fertility receive more fertilizer and a higher irrigation start threshold. This realizes the on-demand allocation of water and fertilizer resources in space, reducing the redundant investment of resources in low soil fertility areas and the risk of resource shortage in high soil fertility areas. Attached Figure Description

[0015] The invention will now be further described with reference to the accompanying drawings.

[0016] Figure 1 This is a flowchart of the system of the present invention; Figure 2 This is a flowchart of the process for generating an intrinsic conductivity distribution map of the entire field in this invention; Figure 3 This is a flowchart of the process for generating the corrected production distribution map in this invention. Detailed Implementation

[0017] The technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

[0018] Please see Figure 1 As shown, the present invention is a smart agricultural control system based on big data, comprising: The data acquisition module is used to collect time-series data of soil electrical conductivity for multiple consecutive growing seasons, historical yield data, elevation data, meteorological data for each grid cell, as well as rice variety attributes and cultivation method attributes corresponding to each growing season. The electrical conductivity denoising module separates the soil electrical conductivity time series data into multiple fluctuation components, removes the high-frequency disturbance components and retains the low-frequency trend components, and generates an intrinsic electrical conductivity distribution map of the entire field by spatial interpolation of the low-frequency trend components. The yield correction module uses the historical yield data of each grid cell as the dependent variable, rice variety attributes and cultivation method attributes as processing variables, and intrinsic conductivity distribution map, elevation data and meteorological data as covariates. It separates the variety effect and cultivation effect through causal inference, and subtracts the variety effect and cultivation effect from the historical yield data to obtain the corrected yield distribution map. The management area division module overlays the intrinsic conductivity distribution map, the corrected yield distribution map, and the elevation data to form a multidimensional feature for each grid. It then compresses the multidimensional feature into a low-dimensional embedding space, simultaneously optimizes clustering allocation and feature reconstruction within the low-dimensional embedding space, and applies embedding similarity constraints to adjacent grids in the space to generate a fine-grained management area boundary. The prescription map generation module determines the unique fertilization amount and irrigation threshold for each management zone based on the median corrected yield and median intrinsic conductivity within each intensive management zone, and generates a variable operation prescription map.

[0019] The data acquisition module is used to collect time-series data of soil electrical conductivity for multiple consecutive growing seasons, historical yield data for each grid cell, elevation data, meteorological data, and rice variety attributes and cultivation method attributes corresponding to each growing season, specifically including: This invention uses paddy fields as an example. Within the target paddy field, several grid units are divided according to a 10m x 10m grid spacing. A soil conductivity sensor is buried at the center of each grid unit, with the sensor probe buried at two depths: 10cm and 20cm. Soil conductivity data is continuously collected every 30 minutes from the date of rice transplanting or direct seeding, covering the entire growing season until harvest. This process is repeated for three consecutive growing seasons, and a set of time-series conductivity data including timestamps is obtained for each sampling point in each growing season.

[0020] During the rice harvesting process in each growing season, the combine harvester's built-in yield monitoring system, which integrates a GPS module and an impulse flow sensor, automatically records a yield data point every 2 seconds, along with its corresponding latitude and longitude coordinates. After harvesting, each yield data point is spatially merged into a 10m x 10m grid, and the average of all yield data points falling within the same grid is taken as the historical yield data for that grid unit in that growing season. This process is repeated to obtain historical yield data for each grid unit across three consecutive growing seasons.

[0021] Within the target rice paddy, a quadcopter drone equipped with lidar was used for low-altitude aerial photography at a height of 50 meters and a flight path overlap rate of 70%. After acquiring point cloud data, a digital elevation model was generated, and the elevation value of each grid cell was extracted. An integrated small weather station was set up in the center of the paddy or in an open area. This station recorded air temperature, humidity, precipitation, wind speed, and solar radiation every hour. All meteorological data was uploaded to a local server in real time via a wireless transmission module.

[0022] Before the start of each growing season, the rice variety name planted in each grid unit is recorded by visiting farmers or reviewing field management records, and categorized as either indica or japonica rice. The cultivation method for each grid unit, including direct seeding or transplanting, is also recorded. These rice variety and cultivation method attributes are then associated with the corresponding grid unit and growing season identifier to form an attribute label dataset that spatiotemporally matches historical yield data.

[0023] Please see Figure 2 As shown, in the conductivity denoising module, the soil conductivity time series data is separated into multiple fluctuation components. After removing the high-frequency disturbance components, the low-frequency trend components are retained. The low-frequency trend components are then spatially interpolated to generate the intrinsic conductivity distribution map of the entire field, specifically including: Adaptive time-frequency decomposition is performed on the soil electrical conductivity time series data obtained at each sampling point; all local maxima and local minima in the electrical conductivity time series data of that sampling point are identified, where a local maximum is a data point whose value is greater than the two adjacent values ​​before and after it, and a local minimum is a data point whose value is less than the two adjacent values ​​before and after it.

[0024] Using cubic spline interpolation, all local maxima are connected by a smooth curve to form the upper envelope, and all local minima are connected by another smooth curve to form the lower envelope.

[0025] The specific method of cubic spline interpolation is to use a cubic polynomial between every two adjacent extreme points so that the entire envelope has continuous first and second derivatives at the extreme points.

[0026] Calculate the average value of the upper and lower envelopes at the same time point to obtain the mean curve. Subtract this mean curve point by point from the original conductivity time series data to obtain the first decomposed component. Repeat the above envelope construction and mean curve calculation process for the remaining residual signal after subtracting the mean curve, obtaining a new decomposed component each time, until the remaining residual signal becomes a monotonic curve or its fluctuation amplitude is less than one-thousandth of the standard deviation of the original data. Through the above steps, multiple decomposed components are obtained.

[0027] For each component, count the number of times it crosses zero within a unit of time, i.e., one hour. Specifically, find the time points on the time series curve of the component where all values ​​change from positive to negative or from negative to positive, record each such crossing event, and accumulate the number of times it occurs within one hour.

[0028] The average period of each decomposition component is calculated, and the average period is equal to two divided by the number of zero crossings per unit time. Decomposition components with an average period shorter than two hours are identified as high-frequency disturbance components and removed from the decomposition component set; decomposition components with an average period longer than twenty-four hours are retained as low-frequency trend components characterizing the intrinsic soil fertility. If the average period of a certain decomposition component falls exactly between two hours and twenty-four hours, its retention is determined based on the ratio of its amplitude to the energy of its adjacent components. Specifically, if the peak amplitude of the decomposition component is less than one-tenth of the average amplitude of all retained low-frequency trend components, it is removed; otherwise, it is retained.

[0029] For each sampling point, all retained low-frequency trend components are summed to obtain the low-frequency trend signal for that sampling point.

[0030] Ordinary Kriging interpolation was used to interpolate the low-frequency trend signal values ​​of each sampling point to every grid cell in the entire field. The specific steps of ordinary Kriging interpolation are as follows: First, the spatial distance between all pairs of sampling points is calculated, and a scatter plot of the experimental variogram is drawn with the distance as the x-axis and half the square of the difference between the low-frequency trend signal values ​​of each pair of sampling points as the y-axis. A spherical model is used to fit the experimental variogram. The parameters of the spherical model include the nugget constant, the sill value, and the range. The nugget constant is the intercept of the variogram when the distance is zero, the sill value is the y-axis value when the variogram tends to be stationary, and the range is the distance corresponding to 95% of the sill value. Using the fitted variogram parameters, a system of Kriging equations is established. For each grid cell to be interpolated, the low-frequency trend signal values ​​of the twelve closest sampling points to that grid cell are selected, and the equations are solved to obtain the weighting coefficients for each sampling point. The low-frequency trend signal values ​​of each sampling point are then weighted and summed according to the weighting coefficients to obtain the interpolation result for that grid cell.

[0031] Repeat the Kriging interpolation process described above, performing interpolation calculations on all grid cells within the entire field sequentially, ultimately generating an intrinsic conductivity distribution map covering the entire field. The value of each grid cell in this distribution map is the stable conductivity background value obtained after high-frequency disturbance removal and spatial interpolation, used for subsequent yield data correction and management zone delineation.

[0032] Please see Figure 3 As shown, in the yield correction module, the historical yield data of each grid cell is used as the dependent variable, rice variety attributes and cultivation method attributes are used as processing variables, and intrinsic conductivity distribution map, elevation data, and meteorological data are used as covariates. The variety effect and cultivation effect are separated through causal inference, and the variety effect and cultivation effect are subtracted from the historical yield data to obtain the corrected yield distribution map, specifically including: For each grid cell, combine the intrinsic conductivity distribution, elevation, and precipitation and accumulated temperature from the meteorological values ​​to form a three-dimensional feature point. The meteorological value is the sum of the cumulative precipitation and the daily average temperature above zero degrees Celsius during the growing season, from the date of transplanting or direct seeding to the date of harvest. Calculate the standardized Euclidean distance between the three-dimensional feature point of the current grid cell and the three-dimensional feature points of all other grid cells within the field.

[0033] The standardized Euclidean distance is calculated as follows: First, calculate the standard deviations of the intrinsic conductivity distribution, elevation, cumulative precipitation, and accumulated temperature of all grid cells. Then, divide the difference in intrinsic conductivity between the current grid cell and another grid cell by the standard deviation of intrinsic conductivity, the difference in elevation by the standard deviation of elevation, the difference in cumulative precipitation by the standard deviation of cumulative precipitation, and the difference in accumulated temperature by the standard deviation of accumulated temperature. Finally, sum the squares of these four standardized differences and take the square root to obtain the standardized Euclidean distance between the two cells. Sort the cells by distance from smallest to largest and select the ten nearest neighboring grid cells as the current grid cell.

[0034] The number of grid cells with the same rice variety attribute as the current grid cell (both are indica rice or both are japonica rice) among the ten selected neighboring grid cells is counted. This number is divided by ten, and the resulting ratio is used as the probability value of the current grid cell receiving the current variety attribute processing.

[0035] Repeat the above process to calculate the probability value of each grid cell in the entire field. Sort all grid cells in ascending order of probability value. Set a difference threshold, which is one-tenth of the range of probability values ​​of all grid cells (maximum value minus minimum value). Starting with the grid cell with the lowest probability value, subtract the probability values ​​of adjacent grid cells one by one. If the difference is less than the threshold, these two grid cells are assigned to the same comparison group; if the difference is greater than or equal to the threshold, a new comparison group is started. Finally, all grid cells in the entire field are divided into several comparison groups, where the probability values ​​of each grid cell in each comparison group are similar, meaning they are similarly likely to receive the current variety attribute treatment.

[0036] Within the same comparison group, the average difference in yield was calculated for different variety attributes. Specifically, the following steps were taken: The historical yield data of all grid cells in the comparison group with the rice variety attribute of indica rice were summed, and the result was divided by the number of indica rice grid cells to obtain the average indica rice yield; the historical yield data of all grid cells in the comparison group with the rice variety attribute of japonica rice were summed, and the result was divided by the number of japonica rice grid cells to obtain the average japonica rice yield; the difference between the average indica rice yield and the average japonica rice yield was used as an estimate of the variety effect.

[0037] Within this comparison group, the average difference in yield corresponding to different tillage attributes was calculated: the historical yield data of all grid cells with direct seeding were summed and divided by the number of direct seeding grid cells to obtain the average yield of direct seeding; the historical yield data of all grid cells with transplanting were summed and divided by the number of transplanting grid cells to obtain the average yield of transplanting; the difference between the average yield of direct seeding and the average yield of transplanting was used as the estimate of the tillage effect.

[0038] For each grid cell, obtain the rice variety attribute and cultivation method attribute actually used in the current growing season. If the variety attribute of the grid cell is indica rice, subtract the above variety effect estimate from the historical yield data of the grid cell; if it is japonica rice, add the variety effect estimate (or understand it as subtracting the negative effect value).

[0039] If the tillage method is direct seeding, the estimated tillage effect is subtracted from the historical yield data; if it is transplanting, the estimated tillage effect is added. The yield value adjusted for variety and tillage effects is used as the corrected yield for that grid cell. This process is repeated for all grid cells to obtain the corrected yield for each grid cell in the entire field. All corrected yield data constitute a corrected yield distribution map. This map eliminates yield differences caused by different varieties and tillage methods, reflecting the yield potential determined solely by soil fertility conditions.

[0040] In the management area delineation module, the intrinsic conductivity distribution map, the corrected yield distribution map, and elevation data are overlaid to form multidimensional features for each grid. These multidimensional features are then compressed into a low-dimensional embedding space. Within this low-dimensional embedding space, clustering allocation and feature reconstruction are simultaneously optimized, and embedding similarity constraints are applied to adjacent grids to generate a refined management area boundary. Specifically, this includes: A three-layer nonlinear compression transform is constructed to map the multidimensional features of each grid cell to a low-dimensional embedding point. The multidimensional features consist of three dimensions: intrinsic conductivity distribution, corrected yield distribution, and elevation. The first layer of the compression transform uses an intermediate hidden layer with five nodes. The output value of each node is the weighted sum of the three dimensions of the input features plus a bias, then mapped through a nonlinear transformation function that takes the hyperbolic tangent of the input values. The second layer of the compression transform takes the five node values ​​from the first layer as input and maps them to a second hidden layer with three nodes. Similarly, a weighted sum plus a bias is used, followed by taking the hyperbolic tangent. The third layer of the compression transform maps the three node values ​​of the second hidden layer to a final low-dimensional embedding point. The output of this embedding point does not undergo a nonlinear transformation; it directly takes the algebraic sum of the weighted sum and the bias. The weighting coefficients and biases in the three-layer compression transform are parameters to be optimized, and their initial values ​​are random numbers with a mean of zero and a standard deviation of 0.01.

[0041] In each iteration, the soft assignment probability from each low-dimensional embedding point to each candidate cluster center is first calculated. There are K candidate cluster centers, and the value of K is determined as follows: All low-dimensional embedding points of the grid cells are sorted by numerical value and divided into five equal groups. The variance of each group is calculated, and the number of embedding points corresponding to the group with the smallest variance is taken as the initial value of K. This initial value is then compared with the sixth value, and the smaller value is taken as the final value of K. For the ... The low-dimensional embedding point of each grid cell, denoted by its coordinates. For the first There are 1 candidate cluster centers, denoted by their coordinates as follows: Then the first The embedding point is assigned to the first... Soft assignment probability of each cluster Calculate using the following formula: ;in, Show low-dimensional embedding points With candidate cluster centers The Euclidean distance between them. Soft assignment probability. The value range is between 0 and 1, and for a fixed value... ,all of The sum is 1.

[0042] While calculating the soft assignment probability, each low-dimensional embedding point is inversely reconstructed back to the original 3D multidimensional feature. The reconstruction process employs a three-layer expansion transform that is mirror-symmetric to the compression transform: the first layer takes a low-dimensional embedding point as input, mapping it to a hidden layer of three nodes, with the weighting coefficients being the transpose of the weighting coefficients of the third layer of the compression transform, and the biases initialized separately; the second layer maps the three nodes to five nodes; the third layer maps the five nodes to three nodes, outputting the reconstructed feature, whose dimension is the same as the original multidimensional feature. The reconstruction similarity between the original multidimensional feature and the reconstructed feature is calculated. Let the original multidimensional feature be at the π / 2 ... The value in each dimension is The reconstructed feature values ​​in the same dimension are Weight coefficients for each dimension Taking the reciprocals of the variances of intrinsic conductivity, corrected output, and elevation over the global range, the reconstruction similarity R is calculated using the following formula: ;in, For the first The absolute difference in dimensions, multiplied by the weight. Then sum them up, add 1, and take the reciprocal to get the result. . The closer the value is to 1, the higher the reconstruction fidelity. To maximize... To achieve this, the gradient descent method is used to reverse-correct the weighting coefficients and biases of each layer of the compression transform. In each iteration, the correction step size of the weighting coefficients is 0.01, and the correction step size of the biases is 0.001, so that the reconstruction similarity increases monotonically in each iteration.

[0043] For spatially adjacent grid cells, the Euclidean distance between their low-dimensional embedding points is calculated, and this distance is multiplied by a spatial weighting coefficient and then added to the clustering loss. Spatially adjacent grid cells are defined as those whose absolute differences in row and column indices in a field are equal to 1. For each pair of adjacent grid cells, the low-dimensional embedding point of one of them is selected. and another low-dimensional embedding point Calculate the Euclidean distance between the two. The spatial weight coefficient λ is taken as the ratio of the current iteration number to the total number of iterations (preset to 500), multiplied by 0.1. The spatial constraint loss is obtained by summing the Euclidean distances of all adjacent grid pairs multiplied by λ. The total clustering loss consists of two parts: first, the KL divergence between the soft assignment probabilities of all grid cells and the target distribution; and second, the aforementioned spatial constraint loss. The target distribution is calculated as follows: [The text abruptly ends here, so the translation stops.] Square it and then divide by the cluster. All grid cells within The sum is then normalized for each grid cell.

[0044] The iterative process continues. After each iteration, the change in the soft assignment probability of each grid cell relative to the previous iteration is calculated. Specifically, this is done by summing the absolute values ​​of the probability differences between all cluster centers of all grid cells and dividing by the product of the total number of grid cells and the number of cluster centers. A preset change threshold of 0.001 is used. Iteration stops when the change calculated in two consecutive iterations is less than this threshold. At this point, for each grid cell, the category label corresponding to the cluster center with the highest soft assignment probability is taken as the category to which the grid cell belongs. Grid cells with the same category label and spatially adjacent cells are merged into a connected region, and the boundary of each connected region is the boundary of the fine-grained management area. All fine-grained management area boundary lines are output for subsequent determination of the fertilization amount and irrigation threshold for each management area.

[0045] In the prescription map generation module, based on the median corrected yield and median intrinsic conductivity within each intensive management zone, the unique fertilizer application rate and irrigation threshold for each zone are determined, generating a variable operation prescription map, specifically including: Obtain the boundaries of all fine-management zones, each containing several grid cells. For each fine-management zone, extract the corrected yield values ​​of all grid cells within that zone, sort these values ​​from smallest to largest, and take the value at the middle position as the median corrected yield for that fine-management zone. Similarly, extract the intrinsic conductivity values ​​of all grid cells within that zone, sort them, and take the median value as the median intrinsic conductivity for that fine-management zone. Sort all fine-management zones in the entire field from highest to lowest according to the median corrected yield, and select the fine-management zone with the highest value as the benchmark reference zone.

[0046] For each non-reference area outside the reference area, two ratios are calculated. The first ratio, called the first ratio, is calculated by dividing the median corrected output of the non-reference area by the median corrected output of the reference area. The second ratio, called the second ratio, is calculated by dividing the median intrinsic conductivity of the non-reference area by the median intrinsic conductivity of the reference area. Both ratios are dimensionless positive numbers.

[0047] For each non-benchmark reference area, the first ratio and the second ratio are added together and divided by two to obtain a comprehensive coefficient. This comprehensive coefficient reflects the overall ratio of the non-benchmark reference area relative to the benchmark reference area in terms of soil fertility potential and basic soil electrical conductivity. A standard fertilizer application rate for the benchmark reference area is pre-set. This standard fertilizer application rate is determined by the local agricultural technology department based on the target rice yield and basic soil fertility, recommending a total nitrogen, phosphorus, and potassium application rate per acre. For example, 15 kg of pure nitrogen, 6 kg of phosphorus pentoxide, and 8 kg of potassium oxide per acre. Simultaneously, a benchmark irrigation threshold is pre-set for the benchmark reference area, which is set to -30 kPa of soil water potential. Irrigation is initiated when the soil water potential falls below this value. The comprehensive coefficient of the non-benchmark reference area is multiplied by the standard fertilizer application rate of the benchmark reference area to obtain the fertilizer application rate for the non-benchmark reference area. The benchmark irrigation threshold of the benchmark reference area is then subtracted from the product of the comprehensive coefficient and the benchmark irrigation threshold to obtain the irrigation threshold for the non-benchmark reference area.

[0048] The location boundaries of each intensive management zone are stored as a vector polygon coordinate sequence. The calculated fertilizer application rate (including the application rate of nitrogen, phosphorus, and potassium per acre) and irrigation threshold (expressed as soil water potential in kPa) are used as the operational parameters for that zone. The boundaries of all intensive management zones are associated with their corresponding operational parameters to form a complete variable operation prescription map. This prescription map can be exported in a data format compatible with variable fertilizer applicators and intelligent irrigation controllers to guide differentiated field operations.

[0049] The working principle of this invention is as follows: Time-series data of soil electrical conductivity, historical yield data, elevation data, meteorological data, and rice variety attributes and cultivation method attributes corresponding to each growing season are collected for multiple consecutive growing seasons. The electrical conductivity time-series data of each sampling point is separated into multiple fluctuation components through adaptive time-frequency decomposition. High-frequency disturbance components with an average period shorter than two hours are removed, while low-frequency trend components with an average period longer than twenty-four hours are retained. Ordinary Kriging interpolation is used to generate an intrinsic electrical conductivity distribution map for the entire field. Using the intrinsic electrical conductivity distribution value, elevation value, and meteorological value of each grid cell as covariates, the probability value of each grid cell receiving the current variety attribute processing is calculated, and grid cells with similar probability values ​​are grouped into the same comparison group. Within each comparison group, estimates of variety effect and tillage effect are calculated separately. The corresponding effect values ​​are subtracted from historical yield data to obtain a corrected yield distribution map. The intrinsic conductivity distribution map, corrected yield distribution map, and elevation data are superimposed on the three-dimensional features of each grid. A three-layer nonlinear compression transformation is used to compress the features into a low-dimensional embedding space. Within this low-dimensional embedding space, soft-assignment clustering and feature reconstruction are simultaneously optimized, and embedding similarity constraints are applied to adjacent grids. After iterative convergence, the boundaries of the refined management area are output. The region with the highest corrected yield median and intrinsic conductivity median within each refined management area is used as the benchmark reference area. The comprehensive coefficients of each non-benchmark reference area are calculated, thereby determining the unique fertilizer application rate and irrigation threshold for each management area and generating a variable operation prescription map.

[0050] The foregoing has provided a detailed description of one embodiment of the present invention, but this description is merely a preferred embodiment and should not be construed as limiting the scope of the invention. All equivalent variations and modifications made within the scope of the claims of this invention should still fall within the patent coverage of this invention.

Claims

1. A smart agricultural control system based on big data, characterized in that, include: The data acquisition module is used to collect time-series data of soil electrical conductivity for multiple consecutive growing seasons, historical yield data, elevation data, meteorological data for each grid cell, as well as rice variety attributes and cultivation method attributes corresponding to each growing season. The electrical conductivity denoising module separates the soil electrical conductivity time series data into multiple fluctuation components, removes the high-frequency disturbance components and retains the low-frequency trend components, and generates an intrinsic electrical conductivity distribution map of the entire field by spatial interpolation of the low-frequency trend components. The yield correction module uses historical yield data for each grid cell as the dependent variable, rice variety attributes and cultivation method attributes as processing variables, and intrinsic conductivity distribution map, elevation data, and meteorological data as covariates. It separates the variety effect and cultivation effect through causal inference, and subtracts these effects from the historical yield data to obtain the corrected yield distribution map. Specifically, it includes: Based on the intrinsic conductivity distribution, elevation, and meteorological data of each grid cell, the probability value of each grid cell receiving the current rice variety attribute treatment is calculated, and grid cells with similar probability values ​​are grouped into the same comparison group. Within the same control group, the average difference in yield was calculated for different variety attributes and different tillage attributes, and the average difference was used as the estimate of variety effect and tillage effect. By subtracting the estimated effect values ​​corresponding to the actual variety attributes and farming methods of the current year from the historical yield data of each grid cell, a corrected yield distribution map is obtained. The calculation of the probability value of each grid cell receiving the current rice variety attribute processing specifically includes: The intrinsic conductivity distribution, elevation, and meteorological values ​​of each grid cell are combined to form the three-dimensional feature points of the grid cell; Calculate the standardized Euclidean distance between the three-dimensional feature point and the three-dimensional feature points of all other grid cells in the field, and select a preset number of neighboring grid cells in order of increasing distance; The proportion of neighboring grid cells that share the same rice variety attribute as the current grid cell is used as the probability value for the current grid cell to receive the current variety attribute processing. Sort all grid cells according to their probability values ​​from low to high, and group adjacent grid cells whose probability values ​​differ from a preset threshold into the same comparison group; The management area division module overlays the intrinsic conductivity distribution map, the corrected yield distribution map, and the elevation data to form a multidimensional feature for each grid. It then compresses the multidimensional feature into a low-dimensional embedding space, simultaneously optimizes clustering allocation and feature reconstruction within the low-dimensional embedding space, and applies embedding similarity constraints to adjacent grids in the space to generate a fine-grained management area boundary. The prescription map generation module determines the unique fertilization amount and irrigation threshold for each management zone based on the median corrected yield and median intrinsic conductivity within each intensive management zone, and generates a variable operation prescription map.

2. The smart agricultural control system based on big data according to claim 1, characterized in that, The generation of the intrinsic conductivity distribution map for the entire field specifically includes: Adaptive time-frequency decomposition was performed on the soil electrical conductivity time series data of each sampling point to obtain multiple decomposition components arranged from low to high frequency; Calculate the average period of each decomposition component, identify the decomposition components with an average period of less than two hours as high-frequency disturbance components and remove them, and retain the decomposition components with an average period of more than 24 hours as low-frequency trend components. Using the unbiased interpolation method with optimal spatial autocorrelation, the low-frequency trend components of all sampling points are interpolated to each grid cell of the entire field to generate an intrinsic conductivity distribution map.

3. The smart agricultural control system based on big data according to claim 2, characterized in that, The adaptive time-frequency decomposition of the soil conductivity time-series data for each sampling point specifically includes: Identify all local extrema in the conductivity time series data, connect the local maxima and local minima using cubic spline interpolation to form upper and lower envelopes, and calculate the mean curves of the upper and lower envelopes; The mean curve is subtracted successively from the original conductivity time series data. Each subtraction yields a decomposed component. The envelope and mean calculation process is repeated until the remaining signal is monotonic or has no obvious fluctuations, thus obtaining multiple decomposed components. Count the number of times each decomposed component crosses zero within a unit of time, and arrange the decomposed components in order of the number of zero crossings from least to most, thus obtaining the decomposed components arranged from low to high frequency.

4. The smart agricultural control system based on big data according to claim 1, characterized in that, The generation of the fine-grained management area boundary specifically includes: A three-layer nonlinear compression transformation is constructed to sequentially map the multidimensional features of each grid into intermediate hidden layer representations and the final low-dimensional embedding points; In each iteration, the soft assignment probability from each low-dimensional embedding point to each candidate cluster center is calculated. At the same time, each low-dimensional embedding point is reconstructed back to the original multi-dimensional features, and the similarity between the features before and after reconstruction is used as a constraint to adjust the parameters of the compression transformation. For spatially adjacent grids, calculate the Euclidean distance between their low-dimensional embedding points, and then multiply the Euclidean distance by the spatial weight coefficient and add it to the clustering loss. When the change in soft assignment probability between two consecutive iterations is less than a preset threshold, the boundary of the fine management area to which each grid belongs is output.

5. A smart agricultural control system based on big data according to claim 4, characterized in that, The parameter for adjusting the compression transformation by using the similarity between features before and after reconstruction as a constraint specifically includes: Each low-dimensional embedding point is sequentially subjected to a three-layer expansion transformation that is mirror-symmetric to the compression transformation, gradually restoring the dimension until a reconstructed feature with the same dimension as the original multidimensional feature is obtained. Calculate the absolute value of the difference between the original multidimensional features and the reconstructed features in each dimension, and then take the reciprocal of the weighted sum of the absolute values ​​of the differences in each dimension to obtain the reconstructed similarity. With the goal of maximizing the reconstruction similarity, the projection direction and bias of each layer of the compression transformation are corrected in reverse, so that the reconstruction similarity increases monotonically in each iteration.

6. The smart agricultural control system based on big data according to claim 1, characterized in that, The generated variable job prescription diagram specifically includes: The corrected median yield of all finely managed areas in the field is sorted from high to low, and the management area with the highest value is selected as the benchmark reference area. For each non-reference area, calculate the first ratio of the corrected median output of the fine management area to the corrected median output of the reference area, and the second ratio of the median intrinsic conductivity of the fine management area to the median intrinsic conductivity of the reference area. The arithmetic mean of the first ratio and the second ratio is used as the comprehensive coefficient. The standard fertilizer application rate of the benchmark reference area is multiplied by the comprehensive coefficient to obtain the fertilizer application rate of the intensive management area. The benchmark irrigation threshold of the benchmark reference area is subtracted from the product of the comprehensive coefficient and the benchmark irrigation threshold to obtain the irrigation threshold of the intensive management area. The location boundaries of each fine management zone are associated with the corresponding fertilization amount and irrigation threshold to generate a variable operation prescription map.

Images