Crop planting and fertilizing information processing method and system

Abstract

The application discloses a kind of crop planting fertilization information processing method and system, belong to crop planting information processing technical field.The method includes: obtaining the multi-source crop condition data containing soil parameter time series, crop growth observation image and weather forecast data;Empirical mode decomposition is carried out to soil parameter time series, and high-frequency fluctuation component and low-frequency trend component are reconstructed;High-frequency fluctuation component is input into pulse neural network model, and the pulse sequence feature reflecting the short-term mutation risk of soil parameter is output;Low-frequency trend component is spatiotemporally aligned and fused with crop growth observation image, and spatial structure code is extracted using graph convolution network;Mutual information entropy of pulse sequence feature and spatial structure code is calculated, and key collaborative feature pair is screened;Key collaborative feature pair is input into sequence-to-sequence model based on attention mechanism, and preliminary fertilization recommendation atlas containing fertilization time window and nutrient ratio is generated combined with weather forecast data.

Description

Technical Field

[0001] This invention relates to the field of crop planting information processing technology, specifically to a method and system for processing crop planting and fertilization information. Background Technology

[0002] Fertilization decisions during crop cultivation heavily rely on a comprehensive analysis of soil conditions, crop growth, and future weather conditions. For analyzing soil parameter time series, conventional methods often employ statistical regression or frequency domain filtering to extract features. However, these methods struggle to simultaneously separate the long-term evolution of soil nutrients from short-term fluctuations caused by external environmental disturbances. Monitoring crop growth typically relies on multispectral imagery, using vegetation index analysis to determine nutrient deficiencies. However, this approach easily confuses changes in soil background with the crop's physiological responses, lacking a deep understanding of the spatiotemporal correlation mechanisms between the soil and crop systems.

[0003] Existing technologies for processing multi-source heterogeneous agricultural data typically employ simple data stitching or manual interpretation after independent modeling, which has the following drawbacks. Firstly, short-term non-stationary fluctuations in soil parameters are often filtered out as noise, ignoring implicit information about potential mutation risks such as microbial activity outbreaks or rapid nutrient leaching, leading to delayed fertilization timing. Secondly, the lack of a structured expression of the spatiotemporal correspondence between long-term soil fertility trends and crop canopy spectral responses makes it difficult to accurately quantify the intensity of nutrient stress at specific farmland locations due to spatial diffusion from neighboring areas, resulting in insufficient basis for spatially differentiated fertilization allocation. Furthermore, existing methods, when dealing with high-dimensional dynamic features, rely on limited quantitative assessment methods for the degree of synergistic influence between features, failing to adaptively capture key decision-making factors with strong instantaneous correlations.

[0004] To address the issues of discarding short-term soil mutation signals and the difficulty in quantifying the spatiotemporal synergistic mechanism of the soil-crop system, breakthroughs are needed in two dimensions: nonlinear signal decomposition and pulse-driven risk perception, and spatiotemporal pattern fusion and dynamic synergistic screening of the soil-crop system. Summary of the Invention

[0005] The purpose of this invention is to provide a method and system for processing information on crop planting and fertilization, in order to solve the problems in the prior art where the risk of short-term soil mutation is ignored and there is a lack of soil-crop spatiotemporal collaborative analysis methods, resulting in a lag in the fertilization decision-making time window and insufficient spatial configuration accuracy.

[0006] To achieve the above objectives, this invention provides the following technical solution: This invention provides a method and system for processing crop planting and fertilization information, addressing the problems of insufficient response to short-term changes in soil parameters, inadequate fusion of multi-source heterogeneous data, and low precision in fertilization recommendations in existing fertilization decision-making processes. In the method section, multi-source agricultural data for the target plot is acquired, including soil parameter time series, crop growth observation images, and weather forecast data. The soil parameter time series is processed using empirical mode decomposition (EMD) to obtain multiple intrinsic mode function (IMF) components and a residual component. IMF components with frequencies exceeding a preset threshold are reconstructed as high-frequency fluctuation components, and the remaining IMF components are merged with the residual component into a low-frequency trend component. The high-frequency fluctuation component captures instantaneous changes and short-term oscillations in soil parameters, accurately reflecting sudden disturbances in the field microenvironment; the low-frequency trend component retains the macroscopic evolution of soil nutrients, providing a stable basis for judging the long-term nutrient demand trend of crops. The high-frequency fluctuation component is input into a pre-trained spiking neural network model, activating the leaky integral ignition neurons within the model, and outputting pulse sequence features reflecting the risk of short-term changes in soil parameters. Spiking neural networks, with their pulse time coding mechanism, have extremely high response sensitivity to sparse burst events in high-frequency fluctuation components, significantly improving the timeliness and accuracy of identifying short-term mutation risks in soil.

[0007] Low-frequency trend components are spatiotemporally aligned and fused with crop growth observation images to obtain a fused feature tensor. A graph convolutional network is then used to extract spatial topological relationships from this fused feature tensor, yielding a spatial structure code. By strictly aligning the soil low-frequency trend sequence with crop growth images in both time and space, the graph convolutional network can fully learn the intrinsic spatial dependence between adjacent plots on crop growth and soil nutrient supply and demand. This allows the spatial structure code to encompass both the real-time crop growth status reflected by vegetation indices and the potential supply capacity represented by soil nutrient trends. The mutual information entropy between pulse sequence features and the spatial structure code is calculated. Based on the dynamic change rate of the mutual information entropy, key synergistic feature pairs that significantly influence fertilization decisions are selected. Mutual information entropy can unbiasedly measure the strength of the nonlinear statistical correlation between short-term mutation risk signals and spatial growth status signals, automatically filtering redundant and weakly correlated features, retaining only feature combinations truly sensitive to nutrient regulation, thereby reducing the computational complexity of subsequent models and improving recommendation accuracy. Key collaborative features are input into an attention-based sequence-to-sequence model, which, combined with meteorological forecast data, generates a preliminary fertilization recommendation map. This map includes the fertilization time window and nutrient ratio for each grid cell. The attention mechanism allows the model to autonomously focus on the interaction between key features influencing fertilization timing and future changes in meteorological conditions. The sequence-to-sequence structure can generate a complete temporal fertilization plan for all grid cells at once, ensuring the consistency of time windows and the coordination of spatial distribution among the cells in the recommendation map.

[0008] As a preferred technical solution of the present invention, the empirical mode decomposition processing of soil parameter time series includes: extracting all local maxima and local minima, constructing upper and lower envelopes using cubic spline interpolation, and calculating candidate components with the mean of the envelope; repeatedly screening the mean of the envelope until the number of local extrema and the number of zero-crossing points of the candidate component differ by no more than one and the mean of the envelope is zero, separating the candidate component into an intrinsic mode function component; subtracting this component from the original sequence to obtain the residual sequence, and repeating the above operation until the residual sequence is monotonic or the number of extrema is lower than a preset value; calculating the zero-crossing density of each intrinsic mode function component to obtain the instantaneous frequency, summing the components with instantaneous frequencies greater than a preset threshold as high-frequency fluctuation components, and using the sum of the remaining components and the residual component as low-frequency trend components. By combining extreme point interpolation with zero-crossing density, this scheme achieves adaptive multi-scale decomposition of soil parameter sequences. Frequency segmentation does not require preset basis functions, can truly reflect the nonlinearity and non-stationarity of the data itself, and avoids subjective errors introduced by manually setting the decomposition level.

[0009] When processing high-frequency fluctuation components, the spiking neural network loads these components as input current signals into the input layer neurons, with each neuron corresponding to a time step value. Leakage integral ignition neurons accumulate membrane potential according to the membrane potential dynamic equation, firing a pulse and resetting the membrane potential to resting potential when the membrane potential exceeds the ignition threshold. The pulse firing time and number of times for each neuron within a preset time window are recorded to form the original pulse sequence. This sequence is then aligned in the time dimension using a sliding window alignment technique and normalized frequency encoding to generate pulse sequence features. This scheme utilizes the leakage integral characteristic of neurons to naturally decay non-critical perturbations, firing pulses only for strong mutation events exceeding the threshold. It is suitable for filtering weak noise in soil sensor data while retaining mutation signals with actual physiological risks, highlighting critical nodes requiring urgent intervention in fertilization decisions.

[0010] When fusion of low-frequency trend components with crop growth observation images in a spatiotemporal alignment, the crop growth observation images are first subjected to multispectral band extraction to obtain the vegetation index time series. Bilinear interpolation is used to synchronize the low-frequency trend components with the vegetation index series in time. The time-synchronized low-frequency series is mapped to the corresponding pixel regions of the image according to the spatial grid position. Spatial weighted averaging is performed on each pixel region to generate a pixel-by-pixel low-frequency feature map. The low-frequency feature map is then concatenated with the current frame of the vegetation index along the channel dimension to obtain the fused feature tensor. A spatial topology graph is established with each grid cell as a graph node. The edge connection weights are determined by the Euclidean distance between grid cells and the soil texture similarity. A graph convolutional network performs convolution operations on the topology graph in the spectral domain to aggregate the features of neighboring nodes to generate spatial structure encoding. The combination of bilinear interpolation and spatial weighted averaging ensures a one-to-one correspondence between soil point source sampling data and area source remote sensing images in terms of time and space, avoiding feature confusion caused by geometric misalignment. The edge weight design that integrates Euclidean distance and soil texture enables the graph convolutional network to simultaneously take into account the spatial nearest neighbor effect and the homogeneity of soil texture, thereby extracting spatial feature representations that are closer to the actual topological relationships of crop root growth and nutrient transport.

[0011] The process of calculating mutual information entropy and screening key synergistic feature pairs includes: dividing the pulse sequence features and spatial structure encoding into multiple segments according to the same time window; for each time window, calculating the joint probability distribution matrix of the pulse sequence segment and the spatial encoding segment and their respective marginal probability distributions, and then calculating the mutual information entropy value of that window; traversing all windows to obtain the mutual information entropy time series; performing a difference operation on the mutual information entropy time series to obtain the dynamic change rate of mutual information entropy; and identifying feature pairs whose change rate exceeds a preset change threshold as key synergistic feature pairs. The preset change threshold is dynamically generated by multiplying the standard deviation of the mutual information entropy time series by a preset coefficient. This scheme utilizes the statistical distribution characteristics of mutual information entropy itself to form an adaptive threshold, overcoming the limitation of fixed thresholds being difficult to adapt to different plots and crop varieties, and ensuring that the screened features always have statistically significant synergistic changes.

[0012] When generating a preliminary fertilization recommendation map using a sequence-to-sequence model, key collaborative feature pairs are encoded into a key-value pair vector sequence, and weather forecast data is encoded into a query vector sequence according to time steps. Both are input into the encoder. The encoder uses a multi-head self-attention mechanism to calculate the dependency weights between different positions in the key-value pair vector sequence and outputs an encoded feature sequence. The decoder receives the encoded feature sequence and uses a cross-attention mechanism to allow the query vector sequence and the encoded feature sequence to interact, decoding the fertilization operation vector for each grid cell step by step. This operation vector is then mapped to the start and end times of the fertilization time window and the NPK ratio, forming the preliminary fertilization recommendation map. This design captures long-distance dependencies between different soil and crop characteristics through multi-head self-attention and dynamically regulates fertilization operations based on weather forecast conditions through cross-attention, thereby achieving joint decision-making on fertilization time and nutrient ratio, significantly reducing the risk of inappropriate fertilization timing or formula mismatch.

[0013] In a further embodiment of the present invention, after acquiring multi-source agricultural data, the data is preprocessed. Missing values ​​are detected in the soil parameter time series, locating the timestamps and soil parameter types of missing data. For each missing location, an autoregressive moving average model is constructed using a predetermined number of normal data points before and after it for predictive imputation. The model order is automatically optimized using the Akaike Information Criterion. Simultaneously, crop growth observation images undergo radiometric calibration and atmospheric correction, converting the original digital quantization values ​​into surface reflectance values ​​to obtain preprocessed images. Meteorological forecast data is spatially downscaled, interpolating its grid resolution to match that of the crop growth observation images. These preprocessing steps improve the integrity and consistency of the multi-source data. The autoregressive moving average model accurately restores the continuity of soil parameters over time, avoiding interference from missing values ​​in subsequent decomposition and modeling. Radiometric calibration and atmospheric correction eliminate sensor response differences and atmospheric attenuation effects, ensuring the physical interpretability of vegetation index calculations. Downscaling ensures precise matching between meteorological grids and field grids, enhancing the spatial specificity of meteorological forecast elements in fertilization decisions.

[0014] After generating the preliminary fertilization recommendation map, a reinforcement learning feedback correction mechanism can be introduced. A historical fertilization record dataset is acquired, including historical fertilization time, location, amount, and corresponding historical crop yield. The recommended fertilization time window and nutrient ratio for each grid cell are extracted from the preliminary fertilization recommendation map and compared with historical records to calculate the difference vector between the recommended and historical values. This difference vector is input into a deep Q-network of the reinforcement learning agent, which outputs adjustment actions for each grid cell, using the difference vector as the state. These actions include increasing, decreasing, or maintaining the recommendation. After executing the adjustment actions, a crop growth model is used to simulate and predict the adjusted crop yield, and the simulated yield is used as a reward signal to update the deep Q-network parameters. The preliminary fertilization recommendation map is corrected based on the optimal adjustment action output by the updated network, generating the final fertilization operation instruction. This feedback loop allows the fertilization scheme to continuously learn from historical high-yield experiences, the crop growth model simulation ensures the safety of the trial-and-error process, and the deep Q-network optimizes fertilization parameters while seeking to maximize yield, giving the final scheme the dual advantages of current multi-source data-driven and historical data-validated approaches.

[0015] As another optimization method, after generating a preliminary fertilization recommendation map, each grid cell is clustered and grouped according to the start time of the fertilization time window to form fertilization batch groups. For each group, the average nutrient ratio and the total fertilization amount based on the total area are calculated. In the time sequence of the batch groups, the average nutrient ratio, total fertilization amount, and grid cell location boundaries are output as fertilization operation task sheets and sent to the task scheduler of the unmanned fertilization equipment. The task scheduler then plans the travel path and fertilization execution parameters. Through batch-based intensive scheduling, the equipment's round-trip refueling and idling consumption are significantly reduced, improving operational efficiency and fertilization uniformity.

[0016] In addition, an automatic verification scheme is provided for the rationality of fertilizer ratios. The nitrogen, phosphorus, and potassium ratios of each grid cell in the preliminary fertilizer recommendation map are compared element-by-element with the standard crop nutrient requirement curve, calculating the deviation of nitrogen, phosphorus, and potassium. When the deviation of any element exceeds a preset deviation threshold, the target grid cell is located, and the supplementary fertilizer amount is calculated based on the current soil nutrient content and the standard expected requirement. This supplementary fertilizer amount is then added to the original nutrient ratio to generate a corrected fertilizer recommendation map, and the corrected grid cell is output as a differentiated fertilization instruction. This verification step can promptly correct extreme ratio errors caused by model bias or data anomalies, ensuring that the final fertilizer formula meets the physiological needs of the crop's growth and development stage, and preventing fertilizer waste and yield reduction due to excess or deficiency of a single element.

[0017] This invention also provides a crop planting and fertilization information processing system, including a memory, a processor, and a computer program stored in the memory and running on the processor. When the processor executes the computer program, it implements the steps of the aforementioned crop planting and fertilization information processing method. This system integrates complex multi-source agricultural data processing, spiking neural network and graph convolution feature extraction, mutual information filtering, sequence-to-sequence map generation, and multi-level optimization and correction functions. It can dynamically analyze the soil-crop-meteorological coupling relationship at different crop growth stages, automatically generate differentiated fertilization schemes with high spatiotemporal resolution, save fertilizer input while ensuring increased and stable crop yields, and is suitable for large-scale precision agriculture applications.

[0018] The technical effects and advantages provided by the present invention in the above technical solution are as follows: Soil parameter time series are decomposed into multiple intrinsic mode function components using empirical mode decomposition, and then reconstructed into high-frequency fluctuation components and low-frequency trend components based on instantaneous frequency. The high-frequency fluctuation components retain detailed perturbation information of soil parameters during non-stationary changes. This high-frequency fluctuation component is input into a spiking neural network model, utilizing the cumulative firing mechanism of leaky integral ignition neurons to generate pulses only when the perturbation intensity exceeds a threshold. This encoding method transforms weak but frequent fluctuations into sparse pulse events, enabling the extraction of soil microbial activity pulses or nutrient burst signals implicit in noise as time-precise pulse sequence features, achieving highly sensitive capture of short-term soil mutation risks. Compared to traditional filtering methods that treat such fluctuations as noise suppression, the dynamic characteristics of spiking neural networks can trigger responses at the nascent stage of mutation risk, providing an event-driven perception basis for the precise forward shift of the fertilization time window.

[0019] After spatiotemporally aligning and fusing low-frequency trend components with crop growth observation images, a graph convolutional network is used to extract spatial topological relationships from the fused feature tensor, with farmland grid cells as nodes. The graph convolution process aggregates soil trend information and crop spectral information from adjacent plots, forming a spatial structure encoding containing spatial dependencies. After obtaining the pulse sequence features and spatial structure encoding, the mutual information entropy between the two is calculated and its dynamic change rate is analyzed. Using the fluctuation of mutual information entropy as a criterion, synergistic feature pairs with strong instantaneous correlation between short-term soil fluctuations and crop spatial responses are automatically screened. This mutual information synergistic screening method does not rely on a pre-set linear correlation threshold and can capture nonlinear dynamic coupling relationships, effectively preserving features that only jointly influence fertilization decisions within specific spatiotemporal windows. The screened key synergistic features are input into a sequence model based on an attention mechanism to decode and generate fertilization time windows and nutrient ratios. This ensures that the final fertilization recommendation responds to early mutation risks and follows the synergistic constraints of field spatial heterogeneity, improving the accuracy of fertilization in both spatiotemporal dimensions. Attached Figure Description

[0020] To more clearly illustrate the technical solutions in the embodiments of this application or the prior art, the drawings used in the embodiments will be briefly introduced below. Obviously, the drawings described below are only some embodiments recorded in this invention. For those skilled in the art, other drawings can be obtained based on these drawings.

[0021] Figure 1 This is a flowchart of the information processing method for crop planting and fertilization; Figure 2 This is a flowchart of crop growth observation image fusion and spatial structure coding. Figure 3 This is a flowchart of a crop planting and fertilization information processing system. Figure 4 This is a graph showing the relationship between the edge connection weights of grid cells and the Euclidean distance under different soil texture similarities. Figure 5 It is a graph showing the mutual information entropy and dynamic change rate curves of pulse sequence characteristics and spatial structure encoding; Figure 6 It is a scatter plot of the fertilization nutrient ratio adjustment actions. Detailed Implementation

[0022] To make the objectives, technical solutions, and advantages of the embodiments of the present invention clearer, 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, 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.

[0023] See Figure 1 This invention provides a method for processing crop planting and fertilization information, comprising: Multi-source agricultural data for the target plot is acquired, including soil parameter time series, crop growth observation images, and meteorological forecast data. The soil parameter time series records the changes in soil moisture content, nutrient content, and other indicators over time; the crop growth observation images are acquired by UAVs or remote sensing platforms, covering the visible and near-infrared bands; the meteorological forecast data includes predicted values ​​for future periods such as temperature, precipitation, and wind speed.

[0024] The soil parameter time series is subjected to empirical mode decomposition to obtain multiple intrinsic mode function components and a residual component. The components with frequencies exceeding a preset threshold among the multiple intrinsic mode function components are reconstructed as high-frequency fluctuation components, and the remaining component and the residual component are merged into a low-frequency trend component.

[0025] The high-frequency fluctuation component is input into a pre-trained spiking neural network model, which activates the leakage integral ignition neuron of the spiking neural network model and outputs a pulse sequence feature that reflects the short-term mutation risk of soil parameters.

[0026] The low-frequency trend components are spatiotemporally aligned and fused with the crop growth observation images to obtain a fused feature tensor. A graph convolutional network is then used to extract the spatial topological relationships from the fused feature tensor to obtain a spatial structure encoding.

[0027] Calculate the mutual information entropy between the pulse sequence features and the spatial structure encoding, and based on the dynamic change rate of the mutual information entropy, screen out key synergistic feature pairs that influence fertilization decisions.

[0028] The key collaborative feature pairs are input into an attention-based sequence model, and combined with the meteorological forecast data, a preliminary fertilization recommendation map is generated. The preliminary fertilization recommendation map includes the fertilization time window and nutrient ratio for each grid cell.

[0029] Example 1: In practice, the soil parameter time series is subjected to empirical mode decomposition (EMD) to obtain multiple intrinsic mode function (IMF) components and one residual component. Components with frequencies exceeding a preset threshold are reconstructed as high-frequency fluctuation components, and the residual component is merged with the high-frequency trend component. Cubic spline interpolation is performed on all local maxima in the soil parameter time series to construct an upper envelope; cubic spline interpolation is also performed on all local minima to construct a lower envelope. The point-by-point means of the upper and lower envelopes are calculated to obtain the mean curve. The mean curve is subtracted from the soil parameter time series to obtain candidate components. The IMF conditions for the candidate components are tested: the number of local extrema and the number of zero-crossing points of the candidate component differ by no more than one, and at any time point of the candidate component, the mean of the envelope determined by the upper and lower envelopes is zero. If the conditions are not met, the candidate component is used as a new input signal to repeat the screening process of constructing the upper and lower envelopes, calculating the mean curve, and subtracting signals until the conditions are met. At this point, the candidate component is extracted as an intrinsic mode function (EMF) component. The extracted EMF component is subtracted from the soil parameter time series to obtain the residual sequence. The operations of identifying local extrema, constructing the envelope, filtering the mean, and extracting EMF components are repeated on the residual sequence until the residual sequence exhibits a monotonic function shape or the number of local extrema in the residual sequence is less than a preset number (set to 3). When the number of local extrema is less than 3, the envelope cannot be reliably fitted, and the decomposition stops. This results in multiple EMF components arranged from high to low frequency and one residual component.

[0030] The zero-crossing density of each intrinsic mode function (IMF) component is calculated. This density is obtained by counting the number of sign changes of the IMF component between adjacent sampling points and then dividing by the total duration of the IMF component. The instantaneous frequency of the IMF component is then calculated from the zero-crossing density. The instantaneous frequency is equal to the product of the zero-crossing density and half the sampling frequency. A preset threshold of 0.05 Hz is set based on historical soil parameter fluctuation analysis. IMF components with instantaneous frequencies greater than 0.05 Hz are classified as high-frequency components. All high-frequency components are summed to form the high-frequency fluctuation component. All IMF components with instantaneous frequencies not greater than 0.05 Hz are summed with the residual components to obtain the low-frequency trend component.

[0031] High-frequency fluctuation components are used as input current signals and fed into the input layer neurons of a pre-trained spiking neural network model. Each input layer neuron corresponds to the value of the high-frequency fluctuation component at a single time step. The pre-trained spiking neural network model includes an input layer, a hidden layer, and an output layer. The hidden layer consists of multiple leakage integral ignition neurons, and the neurons in the input layer and the leakage integral ignition neurons in the hidden layer are fully connected. Based on the membrane potential dynamic equation of the leakage integral ignition neurons, each leakage integral ignition neuron in the hidden layer performs membrane potential accumulation processing on the connected input current signal. The discrete update process of the membrane potential is expressed as: in, This represents the membrane potential at the current moment, expressed in millivolts. Indicates after a time step The membrane potential after the membrane; This represents the discrete time step, with a value of 1 millisecond. This represents the membrane time constant, reflecting the membrane potential leakage rate. It is set to 20 milliseconds to enable the leakage integral ignition neuron to respond to high-frequency fluctuations with frequencies higher than 0.05 Hz. This represents the resting potential, with a value of -70 mV, which is the stable potential of the leaky integral ignition neuron when there is no input. This represents the input current signal received by the leaky integral ignition neuron at the current moment, determined by the weighted sum of the high-frequency fluctuation components transmitted by the input layer neurons and the connection weights. When the accumulated membrane potential exceeds the ignition threshold... At that time, the leakage integral ignition neuron outputs a pulse signal and resets the membrane potential to the resting potential. Then it enters a refractory period lasting 1 to 2 milliseconds, during which the inner membrane potential remains at a certain value. And it is not driven by the input current signal. Ignition threshold The setting is -50 mV, which allows the potential difference between the resting potential and the ignition threshold to control the neuron's sensitivity to high-frequency fluctuations. The pulse firing time and number of pulses of each hidden layer leakage integral ignition neuron within a preset time window are recorded to form the original pulse sequence. The length of the preset time window is consistent with the duration of the high-frequency fluctuations.

[0032] The pre-training process of the spiking neural network model is as follows: High-frequency fluctuation component samples and their corresponding short-term mutation risk category labels for soil parameters are collected from historical multi-source agricultural data. First, unsupervised pre-training is performed on the connection weights from the input layer to the hidden layer using the pulse temporal dependence plasticity rule. In each pulse event, the weights are adjusted according to the time difference between preceding and following synaptic pulses, enabling the hidden layer neurons to learn typical pulse temporal patterns in the high-frequency fluctuation components. Then, the input-to-hidden layer weights are fixed, and supervised learning is performed on the weights from the hidden layer to the output layer using the risk category labels. A surrogate gradient method is used to map the pulse firing rate of the output layer neurons to the risk category probability, and the output layer weights are updated through backpropagation using cross-entropy loss. After training, all network parameters are fixed and used to generate pulse sequence features during the inference phase.

[0033] After obtaining the original pulse sequence, a sliding window alignment technique was used to align the original pulse sequence in the time dimension: the sliding window length was set to 50 time steps, and the sliding step size was 25 time steps. The original pulse sequence was divided into multiple window segments. Within each window segment, the number of pulse firings of all leaky integral ignition neurons was counted, and the segments were aligned and arranged according to the window start time to eliminate the temporal deviation caused by the response delay of different neurons, resulting in an aligned pulse sequence. The aligned pulse sequence was then normalized for frequency encoding. The number of pulse firings of each leaky integral ignition neuron within each window segment was divided by the product of the sliding window length and the total number of hidden layer neurons, and mapped to the [0,1] interval to generate pulse sequence features.

[0034] Example 2: In specific implementation, please refer to Figure 2The low-frequency trend components are spatiotemporally aligned and fused with crop growth observation images to obtain a fused feature tensor. Multispectral band extraction processing is then performed on the crop growth observation images: the digital quantization values ​​of each pixel in the red and near-infrared bands are read, and the normalized vegetation index (NDVI) is calculated using the red and near-infrared reflectance. The NDVI is calculated according to the formula NDVI = (NIR - RED) / (NIR + RED), where NIR represents near-infrared reflectance and RED represents red reflectance. The NDVI is calculated frame-by-frame for all time-sampled frames covered by the crop growth observation images to obtain a vegetation index time series. Each frame in the vegetation index time series is a two-dimensional matrix, with the number of rows and columns consistent with the number of spatial pixel rows and columns in the crop growth observation images. The time sampling points of the low-frequency trend component are aligned with the time sampling points of the vegetation index time series using bilinear interpolation: The sampling timestamp sets of the low-frequency trend component and the vegetation index time series are obtained. When the time difference between a sampling timestamp of the low-frequency trend component and two adjacent sampling timestamps in the vegetation index time series is less than a preset synchronization tolerance, linear interpolation is performed on the values ​​at the corresponding pixel positions of two adjacent frames in the vegetation index time series to generate a vegetation index value aligned with the sampling timestamps of the low-frequency trend component. Similarly, when the time difference between a sampling timestamp of the vegetation index time series and two adjacent sampling timestamps in the low-frequency trend component is less than a preset synchronization tolerance, linear interpolation is performed on two adjacent values ​​in the low-frequency trend component to generate a low-frequency trend component value aligned with the sampling timestamps of the vegetation index time series. The preset synchronization tolerance is set to half the sampling interval. After bilinear interpolation alignment, a time-synchronized low-frequency sequence and a vegetation sequence are obtained, both with the same number of time sampling points.

[0035] The time-synchronized low-frequency sequences are mapped to corresponding pixel regions of the crop growth observation image according to spatial grid locations: the target plot is pre-divided into regular grids, each grid cell has a geographic coordinate range, and each pixel in the crop growth observation image has geographic coordinates. The geographic coordinate range of the grid cell is spatially mapped to the geographic coordinates of the pixel to determine the set of pixels covered by each grid cell. A spatially weighted average is performed on the low-frequency sequence values ​​within each pixel region to generate a pixel-by-pixel low-frequency feature map: Assuming a pixel region contains m pixels, the set of low-frequency sequence values ​​corresponding to that pixel region at a certain sampling point is { , ,..., The spatially weighted average is calculated by using the reciprocal of the Euclidean distance from the pixel to the center of the grid cell as the weight. The low-frequency sequence values ​​in the set are summed according to their weights and then divided by the total weights. The resulting spatially weighted average is assigned to the low-frequency feature map value of that grid cell at the current time sampling point. This spatially weighted average is performed grid-by-grid cell to obtain a pixel-by-pixel low-frequency feature map containing multiple time frames. The pixel-by-pixel low-frequency feature map is then concatenated with the current frame of the vegetation index time series in the channel dimension: Let the size of the pixel-by-pixel low-frequency feature map be H rows and W columns, and the size of the current frame of the vegetation index time series be H rows and W columns. The pixel-by-pixel low-frequency feature map is considered as one channel, and the current frame of the vegetation index time series is considered as another channel. These are stacked in the channel dimension to form a fused feature tensor of size 2×H×W, where H represents the number of grid cell rows and W represents the number of grid cell columns.

[0036] Each grid cell of the target plot is used as a graph node, and the number of graph nodes equals the total number of grid cells. Each graph node is assigned a node index number. Edge connections are constructed based on the spatial adjacency of adjacent grid cells: if two grid cells share a boundary or a vertex, it is determined that the two grid cells have a spatial adjacency relationship, and an edge connection is established between the corresponding two graph nodes in the graph structure; the edge connection weight of adjacent grid cells in the spatial topology graph is calculated based on the Euclidean distance between grid cells and the soil texture similarity. The edge connection weight calculation method is as follows: in, Represents graph nodes With graph nodes Edge connection weights between them; Represents grid cells With grid cells The Euclidean distance between the center points, in meters. The distance is obtained by calculating the geographic coordinates of the center point of the grid cell using the planar distance formula. This represents the distance attenuation parameter. The standard deviation of the Euclidean distance between the center points of all grid cells is set so that the distance difference can be reflected in the weight changes; This represents the distance term weight adjustment factor. The value is 0.6; This represents the weight adjustment factor for the soil texture similarity term. The value is 0.4. and The value of makes the contribution weight of distance factor in edge connection weight slightly higher than the contribution weight of soil texture similarity factor; Represents grid cells With grid cells Soil texture similarity, The cosine similarity between the soil texture vectors of two grid cells is obtained by calculating the similarity between their soil texture vectors. The soil texture vector consists of three components: sand content, silt content, and clay content. The cosine similarity value ranges from 0 to 1. For grid cell pairs that do not have spatial adjacency, the edge connection weight is set to 0, indicating that there is no edge connection between the graph nodes.

[0037] A graph convolutional network (GCNN) is used to perform spectral domain graph convolution operations on the fused feature tensor, aggregating feature information from neighboring graph nodes and outputting a spatial structure encoding. The GCNN consists of two graph convolutional layers. The first layer receives the fused feature tensor as input (input dimension is 2, the number of channels in the fused feature tensor), and outputs a 64-dimensional feature layer. The second layer receives a 64-dimensional feature layer as input and outputs a 128-dimensional feature layer. The spectral domain graph convolution operation is performed based on the graph Laplacian matrix, which is obtained by subtracting the adjacency matrix from the degree matrix of the spatial topological graph. The elements of the adjacency matrix are the edge connection weights. The degree matrix is ​​a diagonal matrix, with diagonal elements representing the sum of the edge connection weights for each graph node. The propagation rule for each layer of the graph convolutional network is as follows: multiply the feature matrix of the previous layer by the learnable weight matrix, then left-multiply by the normalized graph Laplacian matrix, and finally pass through a linear rectified activation function to obtain the feature matrix of the current layer. After two layers of graph convolutional operations, a spatial structure encoding with a size equal to the total number of graph nodes multiplied by 128 dimensions is obtained. The learnable weight matrix of the graph convolutional network is optimized and determined during the training phase using historical fertilization sample data through a backpropagation algorithm, and the trained weight matrix remains unchanged during the inference phase.

[0038] See Figure 4 The horizontal axis in the figure represents the Euclidean distance between different grid cells within the target plot. The unit is meters, ranging from 0 to approximately 1250 meters; the vertical axis represents the edge connection weight between corresponding grid cells. The figure uses three symbols to represent grid cell pairs in different soil texture similarity ranges: dots represent low soil texture similarity (0 to 0.3), squares represent medium soil texture similarity (0.3 to 0.7), and triangles represent high soil texture similarity (0.7 to 1.0). The dashed lines in the figure represent the theoretical edge connection weight curves calculated independently based on the distance term.

[0039] As shown in the figure, the edge connection weights of all soil texture similarity intervals decrease with increasing Euclidean distance, indicating that distance has a negative correlation with edge connection weights. The edge connection weights of grid cell pairs with high soil texture similarity (▲) are generally high, mainly distributed in the range of 0.3 to 0.95, and are generally greater than 0.6 when the distance is close (from 0 to approximately 500 meters), showing that soil texture similarity has a significant positive contribution to edge connection weights. The edge connection weights of grid cell pairs with medium soil texture similarity (■) are distributed between approximately 0.15 and 0.85, generally below the high similarity interval, but still significantly higher than the low similarity interval. The edge connection weights of grid cell pairs with low soil texture similarity (●) are mainly concentrated in the range of 0.05 to 0.4, and tend to approach 0 after the distance exceeds 400 meters, indicating that the edge connection weights are more significantly affected by distance and are generally lower when the similarity is low.

[0040] The dashed theoretical curve only reflects the influence of distance on edge connection weights, which is significantly lower than the discrete data points corresponding to the high and medium soil texture similarity intervals in reality. This demonstrates that soil texture similarity, as a weighting factor in weight calculation, increases the weight of spatially adjacent edge connections. As the Euclidean distance exceeds approximately 800 meters, all edge connection weights tend to approach lower values, indicating that the spatial adjacency between distant grid cells weakens.

[0041] In summary, this figure visually demonstrates the combined influence of distance and soil texture similarity in the edge connection weight calculation method of Example 2 by showing the distribution relationship between edge connection weight and Euclidean distance in different soil texture similarity intervals. It verifies the spatial topology weight construction mechanism in which distance decay and soil texture similarity complement each other.

[0042] Example 3: In practice, the mutual information entropy between the pulse sequence features and the spatial structure encoding is calculated. Based on the dynamic rate of change of the mutual information entropy, key synergistic feature pairs influencing fertilization decisions are selected. The dimensional structure of the pulse sequence features is the number of neurons multiplied by the total number of time steps. The number of neurons corresponds to the number of leaky integral ignition neurons in the hidden layer of the spiking neural network model, and the total number of time steps is the number of time steps for the high-frequency fluctuation components. The spatial structure encoding contains the same number of time slices as the total number of time steps. The dimension of each time slice is the total number of graph nodes multiplied by 128 dimensions, and the total number of graph nodes is equal to the total number of grid cells in the target plot. The pulse sequence features are divided into multiple pulse sequence segments according to time windows. The time window length is set as follows: the total number of time steps of the pulse sequence features is counted, and the integer obtained by dividing the total number of time steps by 10 and rounding down is used as the time window length. The sliding step size is set to half of the time window length and rounded down. The spatial structure encoding is divided into multiple spatial encoding segments according to the same time window length and sliding step size. The pulse sequence segment and the spatial encoding segment under each time window constitute a time-aligned feature pair. All time-aligned feature pairs are arranged in chronological order to form a feature pair sequence.

[0043] For each time-aligned feature pair in the feature pair sequence, calculate the joint probability distribution matrix: Discretize all pulse frequency values ​​contained in the pulse sequence segment of the feature pair into equal-width bins, with a fixed number of bins of 20. The bin width is the difference between the maximum and minimum values ​​in the pulse sequence segment divided by 20. Similarly, discretize all feature values ​​contained in the spatial coding segment of the feature pair into equal-width bins, with a fixed number of bins of 20. Count the frequencies of the discretized pulse sequence segment values ​​falling in the x-th bin and the discretized spatial coding segment values ​​falling in the y-th bin, forming a joint frequency matrix. Divide each element in the joint frequency matrix by the sum of the total frequencies to obtain the joint probability distribution matrix. , where x represents the discrete box index of the pulse sequence segment, and y represents the discrete box index of the spatially coded segment. The marginal probability distributions of the pulse sequence segments are obtained by summing the rows of the joint probability distribution matrix. The marginal probability distribution of the spatially coded segment is obtained by summing the columns. The mutual information entropy value within this time window is calculated according to the following formula: in, This represents the mutual information entropy value within a time window. This represents the joint probability that the discrete bin index of the pulse sequence segment is x and the discrete bin index of the spatially encoded segment is y. This represents the marginal probability of a discrete bin index x in a pulse sequence segment. The marginal probability of a discrete bin index y in the spatially encoded sub-segment is represented by the summation range covering all bin indices. By traversing all time-aligned feature pairs in the feature pair sequence, the mutual information entropy value within each time window is calculated to obtain the mutual information entropy time series.

[0044] The dynamic change rate of mutual information entropy is obtained by performing a difference operation on the mutual information entropy time series: subtracting two adjacent mutual information entropy values ​​in the mutual information entropy time series yields the difference sequence, which is the dynamic change rate sequence of mutual information entropy. A preset change threshold is dynamically generated by multiplying the standard deviation of the mutual information entropy time series by a preset coefficient: the standard deviation of all mutual information entropy values ​​in the mutual information entropy time series is calculated, and the preset coefficient is set to 1.5. The basis for setting the preset coefficient to 1.5 is that feature pairs whose dynamic change rate of mutual information entropy exceeds 1.5 times the standard deviation are identified as significantly changing feature pairs, reducing the interference of random fluctuations on the screening results. Feature pairs corresponding to time windows in the mutual information entropy dynamic change rate sequence whose values ​​exceed the preset change threshold are identified as key collaborative feature pairs.

[0045] Key collaborative feature pairs are encoded into a key-value pair vector sequence. The encoding method is as follows: for each key collaborative feature pair, the frequency values ​​of all neuron pulses at all time steps in the corresponding pulse sequence segment are stretched into a one-dimensional vector, and the dimensional features of all graph nodes at all time steps in the corresponding spatial encoding segment are stretched into a one-dimensional vector. The two one-dimensional vectors are concatenated end to end to obtain a concatenated feature vector. The concatenated feature vector is passed through a key encoder fully connected layer and a value encoder fully connected layer. The output dimension of the key encoder fully connected layer is 256, and the output dimension of the value encoder fully connected layer is 256, resulting in key vectors and value vectors. Each key collaborative feature pair corresponds to a set of key vectors and value vectors. The key vectors of all key collaborative feature pairs are arranged in order to form a key vector sequence, and the value vectors are arranged in order to form a value vector sequence, together forming a key-value pair vector sequence. The length of the key-value pair vector sequence is equal to the total number of key collaborative feature pairs.

[0046] Meteorological forecast data is encoded into a query vector sequence according to time steps: The meteorological forecast data includes hourly raster data of three meteorological elements—temperature, precipitation, and wind speed—for the future forecast period. The meteorological forecast data of all grid cells within the target area are spatially averaged at each time step to obtain the average value vector of meteorological elements at each time step, with a dimension of 3. The average value vector of meteorological elements at each time step is passed through a fully connected layer of a query encoder to output a query vector with a dimension of 256. The query vectors are arranged in order of time steps to form a query vector sequence, and the length of the query vector sequence is equal to the total number of forecast time steps of the meteorological forecast data.

[0047] The attention-based sequence-to-sequence model comprises an encoder and a decoder. The encoder consists of six stacked encoder layers, each containing a multi-head self-attention sublayer and a feedforward sublayer. Both the multi-head self-attention and feedforward sublayers employ residual connections and layer normalization. In the multi-head self-attention sublayer, the key-value pair vector sequence is taken as input. With eight attention heads, each attention head linearly projects the input key vector sequence to obtain a query matrix, then projects the key vector sequence to obtain a key matrix, and finally projects the value vector sequence to obtain a value matrix. The scaled dot product attention output is then calculated. The outputs of the eight attention heads are concatenated and passed through an output projection matrix to obtain the multi-head self-attention sublayer output. The feedforward sublayer consists of two linear transformation layers. The first linear transformation layer maps the 256-dimensional sequence to 1024-dimensional sequence. After passing through a linearly rectified activation function, the second linear transformation layer maps the 1024-dimensional sequence back to 256-dimensional sequence. The encoder calculates the dependency weights between different positions in the key-value pair vector sequence and outputs an encoded feature sequence. The length of the encoded feature sequence is the same as the length of the key-value pair vector sequence, and the dimension of each position is 256.

[0048] The decoder consists of six stacked decoder layers. Each decoder layer contains a masked multi-head self-attention sublayer, a cross-attention sublayer, and a feedforward network sublayer. All sublayers employ residual connections and layer normalization. The masked multi-head self-attention sublayer receives the query vector sequence and uses a masking mechanism to ensure that each time step only focuses on information from the current time step and previous time steps. The cross-attention sublayer uses the output of the masked multi-head self-attention sublayer as the query and the encoded feature sequence output by the encoder as the key and value to calculate the cross-attention output. The structure of the decoder's feedforward network sublayer is the same as that of the encoder's feedforward network sublayer. The decoder decodes time-step by time, outputting a fertilization operation vector at each time step. The fertilization operation vector has a 6-dimensional dimension.

[0049] The fertilization operation vector is mapped to the start and end times of the fertilization time window and the nitrogen, phosphorus, and potassium ratios in the nutrient proportions: the first two dimensions of the fertilization operation vector are processed by the Sigmoid activation function, representing the probabilities of whether the current time step is the start and end time of fertilization, respectively. If the probability of the start time is greater than 0.5, the current time step is determined to be the start time of the fertilization time window; if the probability of the end time is greater than 0.5, the current time step is determined to be the end time of the fertilization time window. The last three dimensions of the fertilization operation vector are processed by the Softmax activation function to obtain the nitrogen, phosphorus, and potassium ratios, with the sum of the three ratios being 1. The sixth dimension of the fertilization operation vector is an idle dimension. The above decoding output and mapping operations are performed on each grid cell of the target plot. The fertilization time windows and nutrient proportions of all grid cells are organized according to the grid cell order to form a preliminary fertilization recommendation map.

[0050] The training steps for the attention-based sequence-to-sequence model are as follows: Collect key collaborative feature pairs from historical fertilization records, weather forecast data samples, and corresponding fertilization time window and nutrient ratio labels. Supervised learning is employed, with a loss function consisting of two parts: a binary cross-entropy loss for the start and end times of fertilization, and a cross-entropy loss for the nitrogen-phosphorus-potassium ratio. The total loss is the weighted sum of these two losses, with weights set to 1.0 and 1.0 respectively. The Adam optimizer is used, with a learning rate of 0.0001, a batch size of 32, and 100 training epochs. The model parameters with the minimum loss on the validation set are selected as the final inference parameters. After training, the weights of all encoder and decoder layers are fixed to generate a preliminary fertilization recommendation map.

[0051] See Figure 5 In the figure, the horizontal axis represents the time window index, and the vertical axis corresponds to the mutual information entropy of the solid curve on the left. The value represents the dynamic rate of change of mutual information entropy, corresponding to the dashed curve on the right side of the ordinate. Mutual information entropy The value ranges from approximately 0.4 to 1.5, exhibiting an overall trend of initially high and then low. Specifically, within the time window index range of 0 to 40, the mutual information entropy remains at a relatively high level of 1.2 to 1.5, indicating a strong information correlation between pulse sequence features and spatial structure encoding. Subsequently, in the range of 40 to 80, the mutual information entropy decreases significantly, dropping to a moderate level of 0.8 to 1.0, reflecting a weakening of the correlation between features. In the range of 80 to 130, the mutual information entropy slightly rebounds but maintains an overall downward trend, remaining within the range of 0.7 to 1.1. Finally, in the range of 130 to 200, the mutual information entropy further decreases to a low level of 0.4 to 0.8, indicating a reduction in collaborative information between the two features.

[0052] The dynamic change rate of mutual information entropy, represented by the dashed line, fluctuates significantly throughout the time period and is symmetrically distributed between -0.3 and 0.35, exhibiting frequent positive and negative jumps. The peak values ​​of the dynamic change rate repeatedly exceed 0.2, indicating a significant change in the correlation between pulse sequence features and spatial structure encoding at the corresponding time window indices. The trough values ​​of the dynamic change rate are below -0.2, indicating a rapid decay in the information interaction between features.

[0053] This figure illustrates the process of screening feature pairs based on mutual information entropy and its dynamic change rate in Example 3. The intervals with higher mutual information entropy values ​​and their corresponding peak intervals of dynamic change rate can serve as candidate time windows for key synergistic feature pairs, as they indicate significant and dynamically changing mutual information between the pulse sequence features and the spatial structure encoding, showing that both have a significant impact on fertilization decisions at this stage. As the time window index increases, the mutual information entropy gradually decreases, indicating a weakening of the synergistic effect of the feature pairs. Dramatic fluctuations in the dynamic change rate suggest abrupt changes in the interaction between features within certain time windows, facilitating the screening of significantly changing key synergistic feature pairs through thresholding.

[0054] Overall, the figure reveals the mutual information entropy and its dynamic changes between the pulse sequence features output by the spiking neural network and the spatial structure encoding extracted by the graph convolutional network, providing an intuitive basis for the dynamic identification and screening of key collaborative feature pairs in Example 3.

[0055] Example 4: In practice, after acquiring multi-source agricultural data for the target plot, missing value detection is performed on the soil parameter time series to locate the timestamp position of the missing data and the corresponding soil parameter type. The soil parameter time series is sampled at equal intervals, with the sampling interval set according to the sensor configuration. Each sampling point in the soil parameter time series is traversed, checking whether the value corresponding to each soil parameter type at each sampling point is marked as null or exceeds a reasonable value range. Soil parameter types include soil moisture content, total nitrogen content, available phosphorus content, available potassium content, and organic matter content. When the value of a soil parameter type at a sampling point is marked as null, the timestamp position and corresponding soil parameter type of that sampling point are recorded. Similarly, when the value of a soil parameter type at a sampling point exceeds a preset reasonable value range, the timestamp position and corresponding soil parameter type of that sampling point are also recorded, and values ​​exceeding the range are considered missing data. The preset reasonable value ranges are set as follows: the preset reasonable value range for soil moisture content is 0.05 to 0.50 g / kg, the preset reasonable value range for total nitrogen content is 0.2 to 5.0 g / kg, the preset reasonable value range for available phosphorus content is 1.0 to 80.0 mg / kg, the preset reasonable value range for available potassium content is 20.0 to 400.0 mg / kg, and the preset reasonable value range for organic matter content is 3.0 to 60.0 g / kg. The above value ranges are determined based on routine observation and statistical data of the physical and chemical properties of field soils.

[0056] For each timestamp of the missing data, an autoregressive moving average model is constructed using a preset number of normal data points before and after the timestamp of the missing data. The preset number is set to 12 normal data points forward and 12 normal data points backward, for a total of 24 normal data points to form the modeling dataset. The selection of 24 data points allows the autoregressive moving average model to capture the short-term temporal dependence of soil parameters while avoiding parameter estimation bias caused by excessively long lags. The order of the autoregressive moving average model is automatically optimized using the Akaike Information Criterion (AIC) on the modeling dataset. A total of 36 combinations are set, with the autoregressive order p ranging from 0 to 5 and the moving average order q ranging from 0 to 5. For each combination of p and q, the modeling dataset is used to fit the autoregressive moving average model using maximum likelihood estimation. The Akaike Information Criterion value is calculated for each combination, and is obtained according to the formula AIC = 2k - 2ln(L), where k represents the total number of model parameters, k equals the sum of p and q, and L represents the log-likelihood function value of the model. The combination of p and q with the smallest Akaike Information Criterion value is selected as the final order of the autoregressive moving average model. Missing data were predicted and imputed using an autoregressive moving average model of a determined order. The prediction direction employed a weighted average of forward and backward predictions: a value was predicted recursively from the 12 normal data points before the missing data time stamp, and a value was predicted recursively from the 12 normal data points after the missing data time stamp. The weights for both forward and backward predictions were set to 0.5. The weighted average of the two predictions was used as the imputed value to fill the corresponding position in the missing data. This modeling, prediction, and imputed process was executed for each located missing data time stamp to obtain a complete soil parameter time series.

[0057] Radiometric calibration correction is performed on crop growth observation images: The radiometric calibration coefficients of the sensor acquiring the crop growth observation images are obtained. These coefficients include gain and offset coefficients, which are read from the sensor's metadata file. The digital quantization values ​​of the original image are converted to surface reflectance values. The conversion relationship is: surface reflectance value equals the product of the digital quantization value and the gain coefficient, plus the offset coefficient. This conversion is performed pixel-by-pixel to obtain the surface reflectance value image. Atmospheric correction is then performed on the surface reflectance values ​​using dark target subtraction: Pixels in vegetation shadow areas or water bodies in the surface reflectance value image are selected as dark target pixels. The surface reflectance values ​​of these dark target pixels in each spectral band are statistically analyzed. The minimum surface reflectance value of each dark target pixel in each spectral band is taken as the atmospheric path radiance estimate. The atmospheric path radiance estimate for the corresponding band is subtracted from the surface reflectance value of each pixel in each band of the surface reflectance value image to obtain the preprocessed crop growth observation image.

[0058] Spatial downscaling is performed on the weather forecast data to interpolate its grid resolution to match that of the crop growth observation imagery. The weather forecast data comes from a numerical weather prediction model, with a typical original grid resolution of 0.1 degrees latitude and longitude, while the crop growth observation imagery typically has a spatial resolution of 0.0001 degrees latitude and longitude. A bilinear interpolation algorithm is used for spatial downscaling. For each target pixel location in the weather forecast data to be interpolated, four neighboring grid center points surrounding the target are found in the original grid. The distances from the target location to these four neighboring grid center points are weighted proportionally in both the longitude and latitude directions, with higher weights assigned to grid center points closer to the target location. The sum of the weather forecast values ​​for the four neighboring grid center points, weighted according to their weights, is the interpolated weather forecast result for the target pixel location. Bilinear interpolation is performed pixel by pixel. The above spatial downscaling processing is applied to the temperature forecast data, precipitation forecast data, and wind speed forecast data respectively to obtain spatially matched meteorological forecast data. The spatial location of each pixel in the spatially matched meteorological forecast data is completely consistent with the spatial location of the corresponding pixel in the crop growth observation image.

[0059] After generating the preliminary fertilization recommendation map, the grid cells in the map are clustered according to the start time of the fertilization time window. All grid cells in the map are traversed, and the start time value of the fertilization time window for each cell is extracted, represented as an integer in hours. Grid cells with the same start time are grouped into the same fertilization batch group, resulting in multiple batch groups, where grid cells within each batch group have the same start time value.

[0060] For each fertilization batch group, the average nutrient ratio is calculated based on the nutrient ratios of each grid cell within that batch group. The nutrient ratios include three values: nitrogen ratio, phosphorus ratio, and potassium ratio. The average nutrient ratio is calculated as follows: sum the nitrogen ratios of all grid cells within the batch group and divide by the total number of grid cells in that batch group; sum the phosphorus ratios of all grid cells within the batch group and divide by the total number of grid cells in that batch group; sum the potassium ratios of all grid cells within the batch group and divide by the total number of grid cells in that batch group. The total fertilization amount is calculated based on the total area of ​​all grid cells within the batch group and the average nutrient ratio. The area of ​​each grid cell is calculated using the spherical area formula based on the geographic coordinate boundary of the grid cell. The total area of ​​the batch group is obtained by summing the areas of all grid cells. The total fertilization amount is calculated as follows: in, This indicates the total amount of fertilizer applied in a fertilization batch group, in kilograms. This indicates the total area of ​​all grid cells within this fertilization batch group, in hectares; This represents the average nitrogen content of this fertilization batch group, dimensionless; This indicates the average phosphorus content of this fertilization batch, dimensionless; This indicates the average potassium content of this fertilization batch, dimensionless; This represents the baseline value for the standard application rate of nitrogen, set at 150 kg per hectare. This value is determined based on a survey of the baseline nitrogen requirements of crops grown on the target plot. This indicates the baseline value for the standard application rate of phosphorus, set at 60 kg per hectare. This value is determined based on a survey of the baseline phosphorus requirements of crops grown in the target plot. This indicates the baseline value for the standard application of potassium, set at 100 kg per hectare. This value is determined based on a survey of the baseline potassium requirements of the crops grown on the target plot.

[0061] Following the chronological order of fertilization batches, the average nutrient ratio, total fertilizer application amount, and corresponding grid cell location boundaries for each batch are output as fertilization task sheets. The chronological order of fertilization batches is arranged from smallest to largest based on the start time of the fertilization time window, with the earliest batch listed first. Each fertilization task sheet includes a task sheet number, execution time, average nitrogen ratio, average phosphorus ratio, average potassium ratio, total fertilizer application amount, and a list of grid cell location boundaries. The grid cell location boundary list contains the latitude and longitude coordinates of the vertices of all grid cells within the fertilization batch group.

[0062] The fertilization task order is sent to the task scheduler of the unmanned fertilization equipment, which is equipped with a positioning module and a variable fertilization execution mechanism. After receiving the fertilization task order, the task scheduler plans the unmanned fertilization equipment's path according to the execution time sequence and grid cell location boundaries in the task order. The path planning covers the area enclosed by the grid cell location boundary list, and a parallel reciprocating path planning algorithm is used to generate a path point sequence. The fertilization execution parameters of the unmanned fertilization equipment are set according to the average nitrogen, phosphorus, and potassium ratios in the task order, setting the opening level of the variable fertilization mechanism, setting the application rate per unit area according to the total fertilizer application amount, and matching the flow control parameters according to the travel speed and width parameters.

[0063] For each grid cell in the preliminary fertilization recommendation map, the nitrogen, phosphorus, and potassium (NPK) ratio values ​​are extracted. These ratios include the nitrogen, phosphorus, and potassium ratios. The standard crop nutrient requirement curves are pre-stored time series of the standard nitrogen, phosphorus, and potassium nutrient requirement ratios corresponding to different crop growth days. The extracted nitrogen ratio is subtracted element-wise from the corresponding standard crop nutrient requirement ratio in the curve, and the absolute value of the difference is used to determine the nitrogen deviation. Similarly, the extracted phosphorus ratio and potassium ratio are subtracted element-wise from the corresponding standard crop nutrient requirement ratio in the curve, and the absolute value of the difference is used to determine the potassium deviation.

[0064] The preset deviation threshold is set to 0.15. The basis for this setting is that when the ratio of any element in the nitrogen, phosphorus, and potassium fertilizer mix deviates from the standard required ratio by more than 0.15, this deviation indicates a significant risk of nutrient imbalance between the fertilizer ratio and the crop's theoretical needs, requiring supplemental fertilization adjustments. The deviations of nitrogen, phosphorus, and potassium are compared to the preset deviation threshold of 0.15. When the deviation exceeds 0.15 for any of these elements, the corresponding grid cell is identified as a target grid cell with a deviation exceeding the preset threshold, and its grid index number is determined.

[0065] The amount of supplemental fertilizer is calculated based on the current soil nutrient content in the soil parameter time series of the target grid cell and the expected nutrient requirements in the standard crop nutrient requirement curve. The current soil nutrient content in the soil parameter time series includes the current total nitrogen content, current available phosphorus content, and current available potassium content, all read from the sampling point closest to the current time in the complete soil parameter time series. The expected nutrient requirements in the standard crop nutrient requirement curve include the expected uptake of nitrogen, phosphorus, and potassium, obtained from the standard crop nutrient requirement curve based on the current crop growth days index, in kilograms per hectare. The amount of nitrogen supplemental fertilizer is calculated as follows: the expected uptake of nitrogen minus the product of the current total soil nitrogen content and the effective conversion coefficient of soil nitrogen. The effective conversion coefficient of soil nitrogen is taken as 0.6, indicating that 60% of the total soil nitrogen is absorbed and utilized by the current crop. When the calculated result is negative, the amount of nitrogen supplemental fertilizer is set to zero. The calculation method for phosphorus supplementation is as follows: The expected phosphorus absorption is subtracted from the product of the current available phosphorus content in the soil and the effective conversion coefficient of phosphorus in the soil. A soil effective conversion coefficient of phosphorus is set to 0.8, indicating that 80% of the available phosphorus in the soil is absorbed and utilized by the current season's crops. When the calculation result is negative, the phosphorus supplementation amount is set to zero. The calculation method for potassium supplementation is as follows: The expected potassium absorption is subtracted from the product of the current available potassium content in the soil and the effective conversion coefficient of potassium in the soil. A soil effective conversion coefficient of potassium is set to 0.7, indicating that 70% of the available potassium in the soil is absorbed and utilized by the current season's crops. When the calculation result is negative, the potassium supplementation amount is set to zero.

[0066] The supplementary fertilization amounts are overlaid onto the nutrient ratios of the target grid cells in the preliminary fertilization recommendation map. The overlay operation is as follows: the nitrogen supplementary fertilization amount is converted into a nitrogen ratio increment by dividing the nitrogen supplementary fertilization amount by the standard application rate of 150 kg / ha; the phosphorus supplementary fertilization amount is converted into a phosphorus ratio increment by dividing the phosphorus supplementary fertilization amount by the standard application rate of 60 kg / ha; and the potassium supplementary fertilization amount is converted into a potassium ratio increment by dividing the potassium supplementary fertilization amount by the standard application rate of 100 kg / ha. The nitrogen ratio increment is added to the original nitrogen ratio of the target grid cell, the phosphorus ratio increment is added to the original phosphorus ratio of the target grid cell, and the potassium ratio increment is added to the original potassium ratio of the target grid cell. The newly obtained nitrogen, phosphorus, and potassium ratio values ​​are then normalized so that the sum of the three is 1. The superimposed and normalized nutrient ratios are written into the corresponding target grid cells in the preliminary fertilization recommendation map to generate the revised fertilization recommendation map. All target grid cells in the revised fertilization recommendation map that have undergone the above-mentioned supplementary fertilization overlay operation are marked as revised grid cells. The grid index number of the revised grid cell, the original NPK ratio, and the revised NPK ratio are compiled and output as differentiated fertilization instructions.

[0067] After generating the preliminary fertilization recommendation map, a historical fertilization record dataset was also acquired. This dataset includes historical fertilization time, location, amount, and corresponding crop yield. Historical fertilization time is stored in date-timestamp format, location in latitude and longitude coordinate format, amount in pure nitrogen, phosphorus, and potassium per hectare format, and yield in grain yield per hectare format. The recommended fertilization time window and recommended nutrient ratio for each grid cell are extracted from the preliminary fertilization recommendation map. The start time of the recommended fertilization time window is used as the recommended fertilization time. The recommended nitrogen, phosphorus, and potassium application rates are obtained by multiplying the nitrogen ratio by 150 kg / ha, the phosphorus ratio by 60 kg / ha, and the potassium ratio by 100 kg / ha, respectively. Compare the historical fertilization time and amount with the corresponding location in the historical fertilization record dataset: Based on the geographic coordinates of the grid cell, search for the historical fertilization record closest to the center point of the grid cell in the historical fertilization record dataset. The distance threshold is half the length of the diagonal of the grid cell. If no record matching the condition is found, skip the grid cell. If a record matching the condition is found, calculate the time difference between the recommended fertilization time and the historical fertilization time in days, calculate the difference between the recommended nitrogen application rate and the nitrogen application rate in the historical fertilization rate, calculate the difference between the recommended phosphorus application rate and the phosphorus application rate in the historical fertilization rate, and calculate the difference between the recommended potassium application rate and the potassium application rate in the historical fertilization rate. Combine the time difference, nitrogen application rate difference, phosphorus application rate difference, and potassium application rate difference into a four-dimensional difference vector.

[0068] The discrepancy vector is input into the reinforcement learning agent, which internally maintains a deep Q-network. The core architecture of the deep Q-network is as follows: the input layer receives a four-dimensional discrepancy vector and has 4 neurons; the hidden layers consist of three fully connected layers: the first fully connected layer has 128 neurons, the second fully connected layer has 128 neurons, and the third fully connected layer has 64 neurons, all using linear rectified functions as activation functions; the output layer is a fully connected layer with 3 neurons, corresponding to the Q-values ​​of the three adjustment actions: increasing fertilization amount, decreasing fertilization amount, and maintaining the recommended fertilization amount. The deep Q-network takes the discrepancy vector as its state input, performs forward propagation calculations on the state, outputs the Q-values ​​corresponding to the three adjustment actions, and selects the adjustment action with the largest Q-value as the output action using an ε-greedy strategy. The ε-greedy strategy is implemented by generating a random number between 0 and 1. If the random number is less than the exploration rate ε, one of the three adjustment actions is randomly selected as the output action. If the random number is not less than the exploration rate ε, the adjustment action with the largest Q value is selected as the output action. The exploration rate ε gradually decreases from 1.0 to 0.05 with each training round.

[0069] After the adjustment action is executed, the crop growth model is used to simulate and predict the adjusted crop yield. The crop growth model takes the fertilization scheme corresponding to the adjustment action, the corresponding soil parameters in the historical fertilization record dataset, and meteorological data as input, runs a complete growing season simulation, and outputs the simulated predicted crop yield. The simulated crop yield is fed back as a reward signal to the deep Q-network. The reward signal is calculated as follows: if the simulated predicted crop yield increases relative to the historical crop yield, the reward value is set as the percentage increase in yield; if the simulated predicted crop yield decreases relative to the historical crop yield, the reward value is set as the negative percentage decrease in yield; if the action of maintaining the recommended fertilization amount is executed, the reward value is additionally increased by a base value of 0.1 to encourage avoiding unnecessary adjustments when the recommended fertilization scheme is reasonable. The target Q-value is calculated using the reward signal and the next time step state through the Bellman equation. The mean squared error between the target Q-value and the output Q-value in the current state is used as the loss function. The network parameters of the deep Q-network are updated using stochastic gradient descent, with a learning rate set to 0.001. After the training of the deep Q network converges, the exploration rate ε is set to 0 during the inference phase. The output action is selected entirely based on the maximum Q value. The initial fertilization recommendation map is corrected based on the optimal adjustment action output by the updated deep Q network. The correction method is to adjust the fertilization nutrient ratio of the corresponding grid cell in the initial fertilization recommendation map by increasing or decreasing the nitrogen, phosphorus and potassium application rates according to the optimal adjustment action, and generate the final fertilization operation instruction.

[0070] Example 5: In specific implementation, please refer to Figure 3 A crop planting and fertilization information processing system includes a memory, a processor, and a computer program stored in the memory and running on the processor. The memory uses a non-volatile storage medium and stores the computer program, multi-source crop data, historical fertilization record datasets, and standard crop nutrient requirement curve data. The processor calls and executes the computer program stored in the memory.

[0071] When the processor executes the computer program, it acquires multi-source agricultural data for the target plot. This multi-source agricultural data includes soil parameter time series, crop growth observation images, and weather forecast data. The soil parameter time series is collected from a soil sensor network deployed within the target plot and stored in the memory via a data transmission module. The crop growth observation images are collected by a multispectral camera mounted on a UAV at preset intervals and transmitted to the memory. The weather forecast data is obtained from a meteorological service interface and stored in the memory.

[0072] The processor performs empirical mode decomposition (EMD) on the soil parameter time series, obtaining multiple intrinsic mode function (IMF) components and one residual component. Components with frequencies exceeding a preset threshold among the IMF components are reconstructed as high-frequency fluctuation components, and the residual component is merged with the residual component into a low-frequency trend component. During the EMD process, the processor performs cubic spline interpolation to construct upper and lower envelopes, calculates the mean curve, and iteratively extracts the IMF components until the residual sequence is monotonic or the number of extreme points is less than 3.

[0073] The processor inputs high-frequency fluctuation components into a pre-trained spiking neural network model, activating the leaky integral ignition neurons of the spiking neural network model, and outputting pulse sequence features reflecting the short-term mutation risk of soil parameters. The network parameters of the spiking neural network model are trained during the construction phase using historical high-frequency fluctuation component samples and corresponding short-term mutation risk category labels of soil parameters. After training, the network parameters are serialized and stored in memory. During the inference phase, the processor loads the network parameters into memory and performs forward propagation.

[0074] The processor performs spatiotemporal alignment and fusion processing on low-frequency trend components and crop growth observation images to obtain a fused feature tensor. The spatiotemporal alignment and fusion processing includes extracting vegetation index time series from the crop growth observation images using multispectral bands, aligning the time sampling points of the low-frequency trend components and vegetation index time series using bilinear interpolation, generating pixel-by-pixel low-frequency feature maps using spatial grid mapping and spatial weighted averaging, and concatenating the pixel-by-pixel low-frequency feature maps with the current frame of the vegetation index time series along the channel dimension. The processor uses a graph convolutional network to extract spatial topological relationships from the fused feature tensor to obtain spatial structure encoding. The graph convolutional network contains two graph convolutional layers. The processor constructs a spatial topological graph using target plot grid cells as graph nodes and adjacent grid cell spatial adjacency relationships as edges. Edge connection weights are calculated based on the Euclidean distance between grid cells and soil texture similarity. Spectral domain graph convolution operations are performed on the graph Laplacian matrix to aggregate neighboring graph node features.

[0075] The processor calculates the mutual information entropy between pulse sequence features and spatial structure encoding. Based on the dynamic change rate of the mutual information entropy, it selects key synergistic feature pairs that influence fertilization decisions. During the calculation, the processor divides the pulse sequence features and spatial structure encoding into segments according to the same time window to form feature pair sequences. It statistically analyzes the joint probability distribution matrix and marginal probability distribution, calculates the mutual information entropy value for each time window to obtain a mutual information entropy time series, performs a difference operation on the mutual information entropy time series to obtain the dynamic change rate of mutual information entropy, and generates a preset change threshold by multiplying the standard deviation of the mutual information entropy time series by a preset coefficient of 1.5. Feature pairs whose dynamic change rate of mutual information entropy exceeds the preset change threshold are selected as key synergistic feature pairs.

[0076] The processor inputs key collaborative feature pairs into an attention-based sequence-to-sequence model, combining this with weather forecast data to generate a preliminary fertilization recommendation map. This map includes the fertilization time window and nutrient ratio for each grid cell. The attention-based sequence-to-sequence model is trained using historical fertilization records and weather forecast data samples during the construction phase. The trained encoder and decoder weights are then stored in memory. The processor encodes key collaborative feature pairs into a key-value pair vector sequence and weather forecast data into a query vector sequence. It generates the encoded feature sequence through the encoder's multi-head self-attention mechanism and decodes the output fertilization operation vector step-by-step through the decoder's cross-attention mechanism. This fertilization operation vector is mapped to the start and end times of the fertilization time window and the nitrogen, phosphorus, and potassium ratio, forming the preliminary fertilization recommendation map.

[0077] When the processor executes a computer program, each specific implementation step can adopt the corresponding processing procedures described in the foregoing embodiments.

[0078] See Figure 6 The horizontal axis of the graph represents the difference between the recommended nitrogen application rate and the historical nitrogen application rate, in kilograms per hectare (kg / hm²), with a value range of approximately -50 to 50. The vertical axis represents the difference between the recommended phosphorus application rate and the historical phosphorus application rate, also in kilograms per hectare, with a value range of approximately -20 to 30. The legend distinguishes three types of fertilization adjustment actions: a solid black triangle represents "increased fertilization," a gray inverted triangle represents "decreased fertilization," and a light gray dot represents "maintaining the recommended fertilization."

[0079] Observing the scatter plot distribution characteristics, it can be seen that the scatter plots of the "increase fertilizer application" category are mainly concentrated in the negative range of the horizontal axis, ranging from approximately -45 to 0, and the vertical axis is mostly distributed between -20 and 0. This indicates that the recommended nitrogen application rate for the corresponding grid unit is lower than the historical nitrogen application rate, and the recommended phosphorus application rate is also lower than the historical phosphorus application rate, so the fertilizer application rate needs to be adjusted upward. The scatter plots of the "decrease fertilizer application" category are mainly concentrated in the region where both the horizontal and vertical axes are positive, with the horizontal axis ranging from approximately 0 to 40 and the vertical axis ranging from approximately 0 to 30. This indicates that the recommended nitrogen and phosphorus application rates for the corresponding grid unit are both higher than the historical application rates, and the fertilizer application rate needs to be reduced in actual adjustment to avoid over-fertilization. The scatter plots of the "maintain recommended fertilizer application rate" category are clustered in the region where both the horizontal and vertical axes are close to zero, with the horizontal axis ranging from approximately -5 to 10 and the vertical axis ranging from approximately -5 to 10. This shows that the recommended fertilizer application rate for the corresponding grid unit is basically consistent with the historical fertilizer application rate, and the original fertilization plan should be maintained.

[0080] This figure reflects the fertilization action adjustment results based on the output of the reinforcement learning agent in Example 5. The difference in nitrogen and phosphorus application rates on the horizontal and vertical axes reflects the difference between the recommended fertilization plan and historical fertilization records. The three symbol forms correspond to the optimal adjustment actions determined by the deep Q-network. The spatial distribution trend of the scatter points clearly distinguishes the grid cell categories that need to increase, decrease, and maintain the fertilization amount, indicating that the system can reasonably generate fertilization adjustment decisions based on historical data and the current recommended plan, assisting in the formation of more accurate final fertilization operation instructions.

[0081] The above description is merely a specific embodiment of this application, but the scope of protection of this application is not limited thereto. Any changes or substitutions that can be easily conceived by those skilled in the art within the scope of the technology disclosed in this application should be included within the scope of protection of this application.

Claims

1. A method for processing information on crop planting and fertilization, characterized in that, include: Acquire multi-source agricultural data for the target plot, including soil parameter time series, crop growth observation images, and meteorological forecast data; The soil parameter time series is subjected to empirical mode decomposition to obtain multiple intrinsic mode function components and a residual component. The components with frequencies exceeding a preset threshold among the multiple intrinsic mode function components are reconstructed as high-frequency fluctuation components, and the residual component and the residual component are merged into a low-frequency trend component. The high-frequency fluctuation component is input into a pre-trained spiking neural network model, which activates the leakage integral ignition neuron of the spiking neural network model and outputs a pulse sequence feature that reflects the short-term mutation risk of soil parameters. The low-frequency trend components are spatiotemporally aligned and fused with the crop growth observation images to obtain a fused feature tensor. A graph convolutional network is then used to extract the spatial topological relationships from the fused feature tensor to obtain a spatial structure encoding. Calculate the mutual information entropy between the pulse sequence features and the spatial structure encoding, and based on the dynamic change rate of the mutual information entropy, screen out key synergistic feature pairs that influence fertilization decisions; The key collaborative feature pairs are input into an attention-based sequence model, and combined with the meteorological forecast data, a preliminary fertilization recommendation map is generated. The preliminary fertilization recommendation map includes the fertilization time window and nutrient ratio for each grid cell.

2. The method for processing crop planting and fertilization information according to claim 1, characterized in that, The process of performing empirical mode decomposition on the soil parameter time series yields multiple intrinsic mode function components and a residual component. Components with frequencies exceeding a preset threshold among the multiple intrinsic mode function components are reconstructed as high-frequency fluctuation components. The remaining component and the residual component are merged into a low-frequency trend component. Cubic spline interpolation is performed on all local maxima and local minima in the soil parameter time series to construct upper and lower envelopes respectively. The mean curves of the upper and lower envelopes are calculated, and the candidate components are obtained by subtracting the mean curves from the soil parameter time series. Repeat the envelope mean screening of the candidate components until the number of local extrema and the number of zero crossings of the candidate components differ by no more than one, and the mean of the envelope is zero. Extract the candidate components that meet the conditions as an intrinsic mode function component. The intrinsic mode function components are subtracted from the soil parameter time series to obtain the residual series. The operation of extracting intrinsic mode function components from the residual series is repeated until the residual series becomes a monotonic function or the number of extreme points is less than a preset number, thus obtaining the plurality of intrinsic mode function components and the residual components. The zero-crossing density of each intrinsic mode function component is calculated, and the instantaneous frequency of each intrinsic mode function component is obtained based on the zero-crossing density. Intrinsic mode function components with instantaneous frequencies greater than the preset threshold are identified as high-frequency components. All high-frequency components are summarized to obtain the high-frequency fluctuation component. All intrinsic mode function components with instantaneous frequencies not greater than the preset threshold are summed with the residual component to obtain the low-frequency trend component.

3. The method for processing crop planting and fertilization information according to claim 1, characterized in that, The high-frequency fluctuation component is input into a pre-trained spiking neural network model, activating the leakage integral ignition neurons of the spiking neural network model, and outputting pulse sequence features reflecting the short-term mutation risk of soil parameters, including: The high-frequency fluctuation component is loaded as an input current signal onto the input layer neurons of the spiking neural network model, and each input layer neuron corresponds to a time step value of the high-frequency fluctuation component. According to the dynamic equation of the membrane potential of the leakage integral ignition neuron, the membrane potential of the leakage integral ignition neuron connected to each input layer neuron is accumulated. When the accumulated membrane potential exceeds the ignition threshold, a pulse signal is output and the membrane potential is reset to the resting potential. Record the pulse firing time and pulse firing number of each leak integral ignition neuron within the preset time window to form the original pulse sequence; The original pulse sequence is aligned in the time dimension using a sliding window alignment technique to obtain an aligned pulse sequence. The aligned pulse sequence is then normalized and frequency encoded to generate the pulse sequence features.

4. The method for processing crop planting and fertilization information according to claim 3, characterized in that, The low-frequency trend components are spatiotemporally aligned and fused with the crop growth observation images to obtain a fused feature tensor. A graph convolutional network is then used to extract spatial topological relationships from the fused feature tensor to obtain spatial structure encoding, including: Multispectral band extraction processing is performed on the crop growth observation images to obtain a vegetation index time series. The time sampling points of the low-frequency trend component are aligned with the time sampling points of the vegetation index time series by bilinear interpolation to obtain a time-synchronized low-frequency sequence and vegetation sequence. The time-synchronized low-frequency sequence is mapped to the corresponding pixel region of the crop growth observation image according to the spatial grid position. The low-frequency sequence values ​​in each pixel region are spatially weighted and averaged to generate a pixel-by-pixel low-frequency feature map. The pixel-by-pixel low-frequency feature map is concatenated with the current frame of the vegetation index time series along the channel dimension to obtain the fused feature tensor. Using each grid cell of the target plot as a graph node, edge connections are constructed based on the spatial adjacency relationship of adjacent grid cells to form a spatial topology graph. The graph convolutional network is used to perform spectral domain graph convolution operation on the fused feature tensor to aggregate the feature information of neighboring nodes and output the spatial structure code.

5. The method for processing crop planting and fertilization information according to claim 4, characterized in that, The edge connection weights of adjacent grid cells in the spatial topology graph are calculated based on the Euclidean distance between the grid cells and the soil texture similarity.

6. The method for processing crop planting and fertilization information according to claim 4, characterized in that, Calculate the mutual information entropy between the pulse sequence features and the spatial structure encoding. Based on the dynamic change rate of the mutual information entropy, select key synergistic feature pairs that influence fertilization decisions, including: The pulse sequence features are divided into multiple pulse sequence segments according to a time window, and the spatial structure encoding is divided into multiple spatial encoding segments according to the same time window, forming a time-aligned feature pair sequence; For each pulse sequence segment and spatially coded segment within a time window, calculate the joint probability distribution matrix, and calculate the marginal probability distributions of the pulse sequence segment and the spatially coded segment respectively. Based on the joint probability distribution matrix and the marginal probability distribution, calculate the mutual information entropy value within the time window, and traverse all time windows to obtain the mutual information entropy time series; The mutual information entropy time series is subjected to a difference operation to obtain the dynamic change rate of mutual information entropy. Feature pairs whose dynamic change rate of mutual information entropy exceeds a preset change threshold are identified as the key collaborative feature pairs.

7. The method for processing crop planting and fertilization information according to claim 6, characterized in that, The preset change threshold is dynamically generated by multiplying the standard deviation of the mutual information entropy time series by a preset coefficient.

8. The method for processing crop planting and fertilization information according to claim 6, characterized in that, The key collaborative feature pairs are input into an attention-based sequence model, and combined with the meteorological forecast data, a preliminary fertilization recommendation map is generated. This preliminary fertilization recommendation map includes the fertilization time window and nutrient ratio for each grid cell, including: The key collaborative feature pairs are encoded into a key-value pair vector sequence, and the weather forecast data is encoded into a query vector sequence according to the time step. Both are input into the encoder of the sequence-to-sequence model based on the attention mechanism. Inside the encoder, a multi-head self-attention mechanism is used to calculate the dependency weights between different positions in the key-value pair vector sequence to generate an encoded feature sequence; In the sequence-to-sequence model decoder based on the attention mechanism, the encoded feature sequence is received, and the query vector sequence is interacted with the encoded feature sequence using a cross-attention mechanism to decode and output the fertilization operation vector of each grid cell step by step. The fertilization operation vector is mapped to the start and end times of the fertilization time window and the nitrogen, phosphorus, and potassium ratio in the fertilization nutrient ratio, and organized into the preliminary fertilization recommendation map according to the grid cells.

9. The method for processing crop planting and fertilization information according to claim 1, characterized in that, The acquisition of multi-source agricultural data for the target plot, including soil parameter time series, crop growth observation images, and weather forecast data, further includes: Missing value detection processing is performed on the soil parameter time series to locate the timestamp position of the missing data and the soil parameter type corresponding to the missing data; For each missing data timestamp position, an autoregressive moving average model is constructed using a preset number of normal data points before and after the missing data timestamp position. The missing data is then predicted and filled using the autoregressive moving average model to obtain a complete soil parameter time series. The crop growth observation images are subjected to radiometric calibration and correction processing. The digital quantization values ​​of the original images are converted into surface reflectance values, and atmospheric correction is performed on the surface reflectance values ​​to obtain preprocessed crop growth observation images. Spatial downscaling is performed on the meteorological forecast data to interpolate the grid resolution of the meteorological forecast data to match the spatial resolution of the crop growth observation image, resulting in spatially matched meteorological forecast data.

10. A crop planting and fertilization information processing system, comprising a memory, a processor, and a computer program stored in the memory and running on the processor, characterized in that, When the processor executes the computer program, it implements the steps of the crop planting and fertilization information processing method as described in any one of claims 1 to 9.

Images