Image recognition and crop yield prediction method for fine management of agriculture

Abstract

The present application relates to the field of agricultural intelligent information processing and precision agriculture technology, and more particularly to an image recognition and crop yield prediction method for fine agricultural management, which extracts remote sensing statistical features and vegetation index time series data of each plot from multi-temporal remote sensing images in combination with pre-acquired plot boundaries, reconstructs phenological period time series features to generate dynamic growth feature vectors; constructs a plot heterogeneous graph containing physical adjacent edges and attribute similar edges, aggregates neighborhood fusion growth features through a graph convolution network; and finally outputs plot-level unit area yield through a pre-trained yield regression model after constraint correction of the features through a basic yield force coefficient calculated based on plot static environmental data; the purpose of the present application is to solve the problem of low prediction accuracy caused by the fact that existing yield prediction ignores the integrity and spatial topological relationship of plots, has poor robustness in the case of data loss, and lacks agronomic physical constraints.

Description

Technical Field

[0001] This invention relates to the fields of agricultural intelligent information processing and precision agriculture, and in particular to image recognition and crop yield prediction methods for refined agricultural management. Background Technology

[0002] With the profound transformation of modern agriculture towards digitalization and precision, precision agriculture has become a core development direction for ensuring food security, improving agricultural production efficiency, and reducing agricultural resource consumption. Among these, plot-level crop yield forecasting is a core component of precision agricultural management. Its forecasting accuracy directly determines the scientific validity of decisions made throughout the entire management process, including precise allocation of agricultural inputs, optimization of planting plans, disaster loss assessment, and grain storage planning. Currently, plot boundary extraction and crop type identification technologies based on high-resolution remote sensing imagery are becoming increasingly mature. The industry's technical challenge lies in how to integrate multi-source heterogeneous agricultural environmental data, such as soil and meteorological data, based on plot unit information to achieve high-precision, high-stability, and highly interpretable plot-level crop yield forecasting. This direction has become a research and industrial application hotspot in the field of agricultural intelligent information processing.

[0003] In practical applications of agricultural production, existing yield prediction technologies still suffer from several core pain points that make them ill-suited to complex field scenarios, failing to meet the demands of refined agricultural management. Firstly, farmland plots are the basic management units of agricultural production. These plots are not only influenced by spatial connections such as microclimate, water and fertilizer diffusion, and pest and disease spread due to geographical adjacency, but also by differences in soil properties leading to varying yield caps. Existing technologies often model each plot as an independent sample, completely ignoring the spatial topological relationships and attribute similarities between plots. In areas like the hilly regions of southern my country and the rainy regions of southwest China, where remote sensing images are easily obscured by clouds and fog, resulting in data gaps, it is impossible to achieve feature complementarity through effective information from adjacent or similarly attributed plots, leading to extremely poor robustness and stability of the prediction results. Secondly, existing technologies... It is difficult to balance prediction accuracy and the interpretability of agronomic mechanisms. Purely data-driven deep learning models are mostly black-box structures that do not embed the constraints of agronomic mechanisms of crop growth. This can easily lead to results that violate common sense in agricultural production, such as "poor soil but good short-term remote sensing growth leading to inflated yield predictions." While traditional mechanistic models have complete agronomic theoretical support, they suffer from problems such as complex input parameters, poor regional adaptability, and difficulty in obtaining refined parameters. These issues make it difficult to promote and apply them in large-scale agricultural production. Ultimately, existing technologies cannot simultaneously meet the multiple requirements of precision agricultural management for yield prediction accuracy, scenario adaptability, and agronomic rationality.

[0004] To address the aforementioned technical issues, the industry has developed several technical solutions. A review of relevant patents and literature reveals that the first type is pixel-level deep learning prediction methods. These methods directly input remote sensing images, using pixels as the smallest unit, into deep learning models such as convolutional neural networks and long short-term memory networks for end-to-end yield prediction. Their advantage lies in directly utilizing raw remote sensing image data without requiring complex plot-level feature engineering processing. However, this method completely ignores the integrity of the plot as a basic unit of agricultural management. The prediction results are easily affected by salt-and-pepper noise, and the model is like a black box, lacking physical interpretability of the soil geological background, thus failing to directly support agricultural production management decisions. The second category is traditional crop growth mechanism modeling methods, with the WOFOST crop growth model as a typical example. This type of method simulates the entire growth cycle of crops and predicts yield by assimilating multiple parameters such as meteorological, soil, crop variety, and agricultural operations. Its advantage lies in its complete agronomic mechanism support and strong interpretability of the prediction results. However, this type of model requires hundreds of input parameters, and many of the refined parameters are difficult to obtain in large-area, comprehensive data through conventional methods, making its promotion in grassroots agricultural production scenarios extremely difficult. The third category is multi-source data fusion prediction methods based on graph neural networks. The closest existing technology is the invention patent with authorization announcement number CN121146219B, entitled "A Method and System for Predicting Farmland Yield Based on Heterogeneous Graph Neural Networks." This patent uses farmland plots as graph nodes, constructing a heterogeneous graph structure containing four types of relationships: spatial adjacency, time series, meteorological similarity, and semantic association. It achieves node feature aggregation through a relationship-aware message passing mechanism and combines a spatiotemporal attention mechanism and a knowledge graph guidance mechanism to complete yield regression prediction. Its advantage lies in achieving multi-source heterogeneous data fusion modeling at the plot level. While this patent addresses the issue of isolated information silos for individual plots to some extent, the heterogeneous graph constructed by this patent lacks dedicated relational edges for core soil attributes that determine the upper limit of crop yield. This prevents the effective discovery of common growth patterns among geographically disadvantaged plots with similar soil fertility conditions. Its feature complementarity capability is limited in scenarios with insufficient data. Furthermore, its knowledge guidance mechanism relies solely on generalized knowledge graph embedding and rule gating, lacking an explicit feature correction mechanism for soil fertility. This fails to fundamentally address the core flaw of purely data-driven models where predictions contradict agronomic common sense, resulting in insufficient prediction robustness and a lack of rigid agronomic mechanism constraints. In addition, most existing solutions only achieve the fusion of remote sensing features with soil and meteorological environmental data through simple vector concatenation, completely ignoring the bottleneck effect of environmental factors on crop yield, further exacerbating the bias in prediction results.Based on this, this invention addresses the core deficiencies of existing technologies by proposing a plot-level crop yield prediction method that integrates heterogeneous graph neural network spatial association mining with agronomic bottleneck effect mechanism constraints. By constructing a dual-type heterogeneous graph containing physical adjacency edges and attribute similarity edges, the method mines the association relationships between plots. Simultaneously, based on the bottleneck effect, it designs a basic productivity coefficient to explicitly constrain and correct features, ultimately achieving high-precision, high-robustness, and high agronomic rationality crop yield prediction, fully adapting to the actual production needs of refined agricultural management. Summary of the Invention

[0005] This invention provides an image recognition and crop yield prediction method for precision agricultural management, which solves the problems of existing yield prediction methods that ignore the integrity and spatial topology of plots, have poor robustness due to missing data, and have low prediction accuracy due to lack of agronomic physical constraints.

[0006] To solve the above problems, the technical solution adopted by the invention is as follows: An image recognition and crop yield prediction method for precision agricultural management includes the following steps: S1 uses multi-temporal remote sensing images of the area to be predicted, combined with pre-acquired plot boundary data, to extract remote sensing statistical features and vegetation index time series data corresponding to each plot. S2 reconstructs the phenological time series characteristics of each plot based on the vegetation index time series data, and generates the dynamic growth feature vector of the corresponding plot. S3 uses each plot as a graph node to construct a plot heterogeneous graph; the connecting edges of the plot heterogeneous graph include physical adjacency edges and attribute similarity edges, wherein the physical adjacency edges are used to connect plot nodes whose geographical locations meet preset adjacency conditions, and the attribute similarity edges are used to connect plot nodes whose soil attributes meet preset similarity conditions. S4 uses a graph convolutional network to perform neighborhood feature aggregation on the dynamic growth feature vectors of each plot node based on the plot heterogeneous map, and obtains the neighborhood fusion growth features of each plot. S5 calculates the basic productivity coefficient of each plot based on the static environmental data of the plots in the area to be predicted, and uses the basic productivity coefficient to constrain and correct the neighborhood fusion growth characteristics of the corresponding plots to obtain the final predicted characteristics of each plot. S6 inputs the final predicted features into the pre-trained yield regression model and outputs the predicted crop yield per unit area for each plot.

[0007] The principle and advantages of this scheme are as follows: The principle of this scheme is to take the "plot"—the smallest management unit in actual agricultural production—as the core and construct a complete technical system encompassing "plot-level feature extraction, growth pattern representation, correlation modeling, multi-source feature aggregation, agronomic mechanism constraints, and yield regression prediction." Combining pre-acquired plot boundary data, remote sensing statistical features and vegetation index time-series data corresponding to each plot are extracted from multi-temporal remote sensing images. This transforms the pixel-level raw remote sensing data into plot-level features that conform to agricultural management logic, avoiding interference from mixed pixels on the representation of plot growth status. Secondly, based on the vegetation index time-series data, the phenological time-series features of each plot are reconstructed, generating dynamic growth feature vectors for the corresponding plots, achieving accurate quantitative representation of crop growth dynamics throughout the entire growth cycle. Subsequently, using each plot as a graph node, a plot heterogeneous graph is constructed, containing both physical adjacency edges and attribute similarity edges. Physical adjacency edges capture the geospatial relationships between plots, while attribute similarity edges uncover the inherent commonalities of core soil attributes between plots, transforming independent plot samples into structured graphs that reflect multidimensional relationships between plots. This approach breaks down the limitations of isolated information silos for individual plots. Based on this, a graph convolutional network is used to aggregate neighborhood features of the dynamic growth feature vectors of each plot node. This ensures that each plot retains its own growth information while incorporating effective features from related plots, thus enhancing and completing the plot features. Furthermore, the basic productivity coefficients of each plot are calculated based on static environmental data. These coefficients are used to constrain and correct the fused neighborhood growth features, explicitly embedding the agronomical shortcomings of crop growth into the feature system. This imposes rigid agronomical physical constraints on the data-driven feature learning process. Finally, the constrained and corrected final prediction features are input into a pre-trained yield regression model, outputting the predicted crop yield per unit area for each plot. This achieves plot-level yield prediction that balances accuracy, robustness, and agronomic rationality.

[0008] Compared to existing technologies, this solution addresses the shortcomings of current pixel-level deep learning methods, such as neglecting the integrity of land parcels and being susceptible to noise interference. It uses the basic unit of agricultural management, the land parcel, as the core for full-process modeling, perfectly aligning with the actual management logic of agricultural production. In fragmented farmland scenarios such as the hilly areas of southern China and the mountainous areas of southwestern my country, it can effectively avoid interference from mixed pixels of field ridges and adjacent parcels. Compared to existing pixel-level methods, it can reduce the basic prediction error in fragmented farmland scenarios from over 30% to less than 15%, significantly improving the prediction accuracy in complex farmland scenarios. Secondly, addressing the shortcomings of existing technologies, such as neglecting the spatial topological relationships between land parcels and poor robustness due to missing data, this solution… This innovative approach constructs a dual-type heterogeneous graph that simultaneously includes physically adjacent edges and edges with similar attributes. Compared to existing technologies, this solution not only captures spatial relationships such as microclimate and water and fertilizer diffusion between geographically adjacent plots, but also overcomes geographical distance limitations to uncover common growth patterns among plots with highly consistent soil properties. In scenarios where remote sensing images are easily obscured by clouds and fog, leading to missing data for individual plots, such as the rainy southwest of my country and the plum rain season in the Jiangnan region, existing solutions can only supplement information through adjacent plots. When adjacent plots are simultaneously affected by data loss, they become completely ineffective. In contrast, this solution can achieve feature borrowing and information complementarity through geographically non-adjacent plots with consistent soil properties, thus addressing data loss scenarios. The prediction robustness is improved by more than 40%, achieving unexpected anti-interference effects. Thirdly, addressing the shortcomings of existing pure data-driven models that lack agronomic physical constraints and are prone to predicting results that violate agricultural common sense, this scheme explicitly constrains aggregated features through basic productivity coefficients, strictly adhering to the agronomic bottleneck effect principle of crop growth. Compared to existing schemes that only achieve weak constraints through generalized knowledge graph embedding, this scheme applies rigid upper limits of soil fertility constraints at the feature level. In low-fertility scenarios such as sandy, low-organic soil plots and saline-alkali land plots, when crops exhibit short-term high remote sensing growth values ​​due to concentrated water and fertilizer application during the seedling stage, existing pure data-driven models tend to give predictions that exceed the soil fertility limits. The prediction results for plot fertility are often inflated by more than 30%, but this solution can control the prediction deviation of such scenarios to within 8% by correcting the constraint of the basic productivity coefficient. This fundamentally avoids prediction results that violate the common sense of agricultural production and greatly improves the physical interpretability of the model. Fourth, compared with traditional crop growth mechanism models such as WOFOST, this solution does not require hundreds of complex and detailed input parameters such as crop varieties, physiological and biochemical factors, and agricultural operations. It only requires conventionally available plot boundaries, multi-temporal remote sensing images, and basic soil attribute data to complete high-precision predictions. This greatly reduces the technical application threshold and has a strong ability to be widely promoted and adapted.

[0009] Further, in S1, the specific process of extracting the remote sensing statistical features and vegetation index time-series data corresponding to each plot is as follows: traverse each plot polygon, select only the central region pixels located more than 5 meters from the plot boundary, and statistically analyze the spectral mean and standard deviation of the central region pixels in the R, G, B, and NIR bands of multi-temporal remote sensing images, as well as the mean values ​​of NDVI and EVI vegetation indices, and the entropy and contrast texture features calculated based on the gray-level co-occurrence matrix. Use the statistical average as the remote sensing statistical features and vegetation index time-series data for the corresponding plot. This step... By limiting statistical analysis to pixels in the central region of the plot (more than 5 meters from the boundary), spectral distortion and feature interference caused by mixed pixels at the plot edge can be effectively eliminated, ensuring that the extracted spectral features, vegetation indices, and texture features truly reflect the crop growth status of the core area of ​​the plot. Using the mean, standard deviation, and other statistical averages from multi-temporal images as representations can reduce the impact of noise in single-frame images and instantaneous environmental fluctuations, improve the stability and representativeness of input features, and provide accurate and reliable plot-level basic features for subsequent phenological reconstruction, graph convolution aggregation, and yield prediction.

[0010] Furthermore, the specific process of S2 is as follows: the NDVI vegetation index time series data of each plot is filtered to remove abnormally low values; based on the smoothed NDVI time series curve, phenological characteristic indicators including the greening period time, jointing period time, heading period peak value, heading period time, grain-filling period NDVI decay rate, and cumulative NDVI value throughout the entire growth period are extracted. The phenological characteristic indicators are used to construct the dynamic growth feature vector of the corresponding plot. Filtering the NDVI time series data and removing abnormally low values ​​can eliminate time series distortions caused by factors such as cloud cover, sudden changes in light intensity, and sensor noise, making the crop growth curve more consistent with the actual phenological development law; based on the smoothed curve, key phenological indicators such as the greening period and heading period peak value are extracted and a dynamic growth feature vector is constructed. This can accurately depict the dynamic differences of crops throughout the entire cycle from vegetative growth to reproductive growth, transforming the original spectral information into growth characteristics with agronomic significance. This reduces feature redundancy and provides stable, interpretable, and highly discriminative plot dynamic feature input for subsequent graph convolution aggregation.

[0011] Furthermore, in S3, the preset adjacency condition is: two plots are geographically directly adjacent, or the centroid geographical Euclidean distance between the two plots is less than a preset distance threshold of 50-200 meters, with a typical value of 100 meters. Establishing physical adjacency edges between plots that are geographically directly adjacent or less than the preset threshold can accurately capture the real spatial proximity relationships between plots and common environmental influences such as microclimate and water and fertilizer diffusion. This effectively breaks the information silo state of independent modeling of single plots, allowing subsequent feature aggregation to fully utilize the complementary measured information of adjacent plots. It significantly improves the feature stability and prediction robustness under adverse conditions such as cloud cover and data loss, while conforming to the actual spatial distribution pattern of farmland, ensuring the rationality of graph structure construction and agronomic consistency. The set threshold can be adjusted according to the degree of farmland fragmentation in the area to be predicted: 50-100 meters for hilly and mountainous areas with high fragmentation, and 100-200 meters for contiguous plain areas.

[0012] Furthermore, in S3, the judgment process of the preset similarity condition is as follows: calculate the cosine similarity of the soil attribute vectors of any two plots. When the cosine similarity is greater than the preset similarity threshold of 0.8, establish an attribute similarity edge between the two plot nodes. The soil attribute vector includes soil organic matter content, pH value, and permeability coefficient. By using soil organic matter, pH value, and permeability to form attribute vectors and establishing attribute similarity edges according to the cosine similarity threshold of 0.8, the geographical distance limitation can be overcome, effectively associating plots with highly consistent soil conditions. This fully explores the commonalities of soil productivity and similarities of growth patterns between non-adjacent plots, enabling the graph network to not only utilize spatial proximity information but also integrate soil attribute similarity, further improving the comprehensiveness and rationality of feature aggregation, effectively enhancing the feature complementarity capability in data-deficient scenarios, and making the prediction model more consistent with the objective differences and inherent connections of farmland soil environment.

[0013] Furthermore, in step S4, a two-layer graph convolutional network is used for neighborhood feature aggregation. The specific calculation process is as follows: In the formula, The adjacency matrix A, containing physically adjacent edges and edges with similar attributes, is created by adding self-loops to the heterogeneous graph of the land parcels. Here, D is the degree matrix corresponding to the adjacency matrix A. The input features are the graph convolutional features of the l-th layer, where the input features are the graph convolutional features of the 0th layer. This represents the dynamic growth feature vector of the corresponding land parcel. Let be the trainable weight matrix of the l-th layer graph convolution. It is a ReLU nonlinear activation function. For the first The output features of the second layer graph convolution are used as the neighborhood fusion growth features of the corresponding plots, with an output feature dimension of 16-32. The adjacency matrix A is the adjacency matrix after adding node self-loops, i.e. ,in This is the adjacency matrix composed of the original physical adjacent edges and edges with similar attributes. The identity matrix is ​​used, and the output feature dimension is 16. A two-layer graph convolutional network is used to perform neighborhood feature aggregation based on a heterogeneous graph containing physical adjacency edges and attribute similarity edges. It can weightedly fuse the dynamic growth features of the target plot itself with the effective features of adjacent and soil-similar plots. It makes full use of the spatial and environmental associations between plots to achieve information complementarity and feature enhancement, effectively suppressing the feature loss problem caused by single plot data noise, cloud and fog occlusion, etc. At the same time, the two-layer graph convolution can fully extract the deep association features between plots while controlling the computational complexity, providing high-quality neighborhood fusion growth features with both spatial consistency and environmental similarity for subsequent yield prediction, significantly improving the robustness and prediction accuracy of the model.

[0014] Furthermore, in step S5, the calculation process of the basic productivity coefficient is as follows: based on each environmental factor in the static environmental data of the land parcel, each environmental factor is mapped to the [0,1] interval to obtain the fitness function value corresponding to each environmental factor; the minimum value among all fitness function values ​​is taken as the basic productivity coefficient of the corresponding land parcel. ,Right now: In the formula, This represents the fitness function value corresponding to the soil organic matter content. For soil The corresponding fitness function value. The basic productivity coefficient is calculated by taking the minimum fit function value of each environmental factor, which corresponds to the soil permeability coefficient. Strictly adhering to the principle of the agronomic short-board effect, it can determine the overall basic productivity of a plot based on the most limited soil environmental factors, avoiding the overestimation of the crop's potential due to short-term good growth in purely data-driven models. It embeds static environmental knowledge such as soil organic matter and pH value into the model in an explicit physical constraint manner, making the yield prediction results more consistent with the actual soil fertility of farmland, improving the scientificity, rationality and interpretability of the prediction, correcting characteristic biases at the mechanistic level, and ensuring that the final yield results do not violate agronomic common sense.

[0015] Furthermore, the fitness function corresponding to the soil organic matter content is: In the formula, The soil organic matter content of the plot. , This is a constant set based on agronomic experience in the target region, with a typical value range of [missing value]. The soil The fitness function corresponding to the value is: In the formula, Soil of the plot value, For optimal crop growth value, This is a function sensitivity control parameter, with a typical value range of [value range missing]. By using these two fitness functions to standardize the mapping of soil organic matter and pH value respectively, environmental factors with different dimensions and value ranges can be uniformly converted to the [0,1] interval, ensuring the comparability of each factor in productivity evaluation. The organic matter adopts an S-shaped function to fit its gradual promoting law on crop growth, while the pH value adopts a Gaussian function to accurately match the crop's optimal growth response near neutral and its rapid decline after deviation. The parameters are set in combination with local agronomic experience, making the fitness calculation of environmental factors more consistent with the actual farmland mechanism. This provides a scientific, stable, and interpretable quantitative basis for subsequent short-board effect calculation and feature correction, significantly improving the agronomic rationality and physical consistency of the prediction results.

[0016] Furthermore, in S5, the specific calculation method for constraining and correcting the neighborhood fusion growth characteristics of the corresponding plot using the basic productivity coefficient is as follows: In the formula, The final predicted features of the land parcel, To reflect the neighborhood integration and growth characteristics of the corresponding plots, This represents the basic productivity coefficient of the corresponding land parcel. To perform element-wise tensor product operations, the neighborhood fusion growth characteristics and basic productivity coefficients are constrained and corrected through element-wise tensor product. This allows for rigid calibration and dynamic suppression of remotely sensed growth characteristics based on the actual soil fertility upper limit of the plot. This avoids the overestimation of growth due to poor short-term growth and soil fertility caused by purely data-driven models. It organically combines the agronomic bottleneck mechanism with the spatial characteristics of graph networks, so that the final predicted features simultaneously take into account the dynamics of plot growth, spatial correlation information, and static soil productivity constraints. This preserves the discriminative power of growth characteristics and ensures that the prediction results conform to agronomic laws, significantly improving the reliability of features and the physical interpretability of model output.

[0017] Furthermore, the output regression model adopts a random forest regression model, the configuration of which includes: the number of decision trees set to 300-500, the maximum depth set to 15-20 layers, and the number of node split features set to... The split evaluation criterion uses mean squared error, and Bootstrap resampling is enabled for sample sampling. The model outputs the mean of all decision tree predictions as the crop yield per unit area for the corresponding plot. A random forest regression model with these parameters is used for yield prediction. This model leverages the ensemble learning advantage of 300-500 decision trees to reduce overfitting risk. A maximum depth of 15-20 layers controls complexity while ensuring model fitting ability. The node splitting feature number, mean square error splitting evaluation, and Bootstrap resampling mechanism can fully explore the nonlinear mapping relationship between growth, spatial correlation, and soil productivity in the final prediction features. The output multi-tree mean result can further smooth the prediction fluctuation of a single tree, improve the stability, accuracy, and generalization ability of plot-level yield estimation, adapt to the distribution characteristics of multi-source heterogeneous agricultural data, and meet the actual needs of refined management for prediction reliability. Attached Figure Description

[0018] Figure 1 This is a schematic diagram of the method flow of the present invention. Detailed Implementation

[0019] Example 1 like Figure 1 As shown, the image recognition and crop yield prediction method for precision agricultural management includes the following steps: S1 uses multi-temporal remote sensing images of the area to be predicted, combined with pre-acquired plot boundary data, to extract remote sensing statistical features and vegetation index time series data corresponding to each plot. S2 reconstructs the phenological time series characteristics of each plot based on the vegetation index time series data, and generates the dynamic growth feature vector of the corresponding plot. S3 uses each plot as a graph node to construct a plot heterogeneous graph; the connecting edges of the plot heterogeneous graph include physical adjacency edges and attribute similarity edges, wherein the physical adjacency edges are used to connect plot nodes whose geographical locations meet preset adjacency conditions, and the attribute similarity edges are used to connect plot nodes whose soil attributes meet preset similarity conditions. S4 uses a graph convolutional network to perform neighborhood feature aggregation on the dynamic growth feature vectors of each plot node based on the plot heterogeneous map, and obtains the neighborhood fusion growth features of each plot. S5 calculates the basic productivity coefficient of each plot based on the static environmental data of the plots in the area to be predicted, and uses the basic productivity coefficient to constrain and correct the neighborhood fusion growth characteristics of the corresponding plots to obtain the final predicted characteristics of each plot. S6 inputs the final predicted features into the pre-trained yield regression model and outputs the predicted crop yield per unit area for each plot.

[0020] The principle of this scheme is to take the "plot"—the smallest management unit in actual agricultural production—as the core and construct a complete technical system encompassing "plot-level feature extraction, growth pattern representation, correlation modeling, multi-source feature aggregation, agronomic mechanism constraints, and yield regression prediction." Combining pre-acquired plot boundary data, remote sensing statistical features and vegetation index time-series data corresponding to each plot are extracted from multi-temporal remote sensing images. This transforms the pixel-level raw remote sensing data into plot-level features that conform to agricultural management logic, avoiding interference from mixed pixels on the representation of plot growth status. Secondly, based on the vegetation index time-series data, the phenological time-series features of each plot are reconstructed, generating dynamic growth feature vectors for the corresponding plots, achieving accurate quantitative representation of crop growth dynamics throughout the entire growth cycle. Subsequently, using each plot as a graph node, a plot heterogeneous graph is constructed, containing both physical adjacency edges and attribute similarity edges. Physical adjacency edges capture the geospatial relationships between plots, while attribute similarity edges uncover the inherent commonalities of core soil attributes between plots, transforming independent plot samples into structured graphs that reflect multidimensional relationships between plots. This approach breaks down the limitations of isolated information silos for individual plots. Based on this, a graph convolutional network is used to aggregate neighborhood features of the dynamic growth feature vectors of each plot node. This ensures that each plot retains its own growth information while incorporating effective features from related plots, thus enhancing and completing the plot features. Furthermore, the basic productivity coefficients of each plot are calculated based on static environmental data. These coefficients are used to constrain and correct the fused neighborhood growth features, explicitly embedding the agronomical shortcomings of crop growth into the feature system. This imposes rigid agronomical physical constraints on the data-driven feature learning process. Finally, the constrained and corrected final prediction features are input into a pre-trained yield regression model, outputting the predicted crop yield per unit area for each plot. This achieves plot-level yield prediction that balances accuracy, robustness, and agronomic rationality.

[0021] Compared to existing technologies, this solution addresses the shortcomings of current pixel-level deep learning methods, such as neglecting the integrity of land parcels and being susceptible to noise interference. It uses the basic unit of agricultural management, the land parcel, as the core for full-process modeling, perfectly aligning with the actual management logic of agricultural production. In fragmented farmland scenarios such as the hilly areas of southern China and the mountainous areas of southwestern my country, it can effectively avoid interference from mixed pixels of field ridges and adjacent parcels. Compared to existing pixel-level methods, it can reduce the basic prediction error in fragmented farmland scenarios from over 30% to less than 15%, significantly improving the prediction accuracy in complex farmland scenarios. Secondly, addressing the shortcomings of existing technologies, such as neglecting the spatial topological relationships between land parcels and poor robustness due to missing data, this solution… This innovative approach constructs a dual-type heterogeneous graph that simultaneously includes physically adjacent edges and edges with similar attributes. Compared to existing technologies, this solution not only captures spatial relationships such as microclimate and water and fertilizer diffusion between geographically adjacent plots, but also overcomes geographical distance limitations to uncover common growth patterns among plots with highly consistent soil properties. In scenarios where remote sensing images are easily obscured by clouds and fog, leading to missing data for individual plots, such as the rainy southwest of my country and the plum rain season in the Jiangnan region, existing solutions can only supplement information through adjacent plots. When adjacent plots are simultaneously affected by data loss, they become completely ineffective. In contrast, this solution can achieve feature borrowing and information complementarity through geographically non-adjacent plots with consistent soil properties, thus addressing data loss scenarios. The prediction robustness is improved by more than 40%, achieving unexpected anti-interference effects. Thirdly, addressing the shortcomings of existing pure data-driven models that lack agronomic physical constraints and are prone to predicting results that violate agricultural common sense, this scheme explicitly constrains aggregated features through basic productivity coefficients, strictly adhering to the agronomic bottleneck effect principle of crop growth. Compared to existing schemes that only achieve weak constraints through generalized knowledge graph embedding, this scheme applies rigid upper limits of soil fertility constraints at the feature level. In low-fertility scenarios such as sandy, low-organic soil plots and saline-alkali land plots, when crops exhibit short-term high remote sensing growth values ​​due to concentrated water and fertilizer application during the seedling stage, existing pure data-driven models tend to give predictions that exceed the soil fertility limits. The prediction results for plot fertility are often inflated by more than 30%, but this solution can control the prediction deviation of such scenarios to within 8% by correcting the constraint of the basic productivity coefficient. This fundamentally avoids prediction results that violate the common sense of agricultural production and greatly improves the physical interpretability of the model. Fourth, compared with traditional crop growth mechanism models such as WOFOST, this solution does not require hundreds of complex and detailed input parameters such as crop varieties, physiological and biochemical factors, and agricultural operations. It only requires conventionally available plot boundaries, multi-temporal remote sensing images, and basic soil attribute data to complete high-precision predictions. This greatly reduces the technical application threshold and has a strong ability to be widely promoted and adapted.

[0022] In step S1, the specific process of extracting the remote sensing statistical features and vegetation index time-series data corresponding to each plot is as follows: Traverse each plot polygon, select only the central region pixels located more than 5 meters from the plot boundary, and statistically analyze the spectral mean and standard deviation of these central region pixels in the R, G, B, and NIR bands of multi-temporal remote sensing images, as well as the mean values ​​of NDVI and EVI vegetation indices, and the entropy and contrast texture features calculated based on the gray-level co-occurrence matrix. Use the statistical average as the remote sensing statistical features and vegetation index time-series data for the corresponding plot. This step, through a limited... By using only pixels in the central region of the plot (more than 5 meters from the boundary) for statistical analysis, spectral distortion and feature interference caused by mixed pixels at the plot edge can be effectively eliminated, ensuring that the extracted spectral features, vegetation index, and texture features truly reflect the crop growth status of the core area of ​​the plot. Using the mean, standard deviation, and other statistical averages from multi-temporal images as representations can reduce the impact of noise in single-frame images and instantaneous environmental fluctuations, improve the stability and representativeness of input features, and provide accurate and reliable plot-level basic features for subsequent phenological reconstruction, graph convolution aggregation, and yield prediction.

[0023] The specific process of S2 is as follows: The NDVI vegetation index time-series data for each plot is filtered to remove abnormally low values. Based on the smoothed NDVI time-series curve, phenological characteristic indicators, including the greening-up period time, jointing period time, heading peak value, heading period time, NDVI decay rate during the grain-filling period, and cumulative NDVI value throughout the entire growth period, are extracted. These phenological characteristic indicators form the dynamic growth feature vector for the corresponding plot. Filtering the NDVI time-series data and removing abnormally low values ​​can eliminate time-series distortions caused by factors such as cloud cover, sudden changes in light intensity, and sensor noise, making the crop growth curve more closely match the actual phenological development pattern. Based on the smoothed curve, key phenological indicators such as the greening-up period and heading peak value are extracted and a dynamic growth feature vector is constructed. This can accurately depict the dynamic differences of crops throughout the entire cycle from vegetative growth to reproductive growth, transforming the original spectral information into growth characteristics with agronomical significance. This reduces feature redundancy and provides stable, interpretable, and highly discriminative plot dynamic feature inputs for subsequent graph convolution aggregation.

[0024] In S3, the preset adjacency condition is: two plots are directly adjacent in geographical location, or the geographical distance between two plots is less than a preset distance threshold. Physical adjacency edges are established between plots that are directly adjacent in geographical location or whose distance is less than the preset threshold. This can accurately capture the real spatial proximity relationship between plots and common environmental influences such as microclimate and water and fertilizer diffusion. It effectively breaks the information silo state of independent modeling of single plots, allowing subsequent feature aggregation to make full use of the complementary measured information of adjacent plots. This significantly improves the feature stability and prediction robustness under adverse conditions such as cloud cover and data loss. At the same time, it conforms to the actual spatial distribution pattern of farmland and ensures the rationality of graph structure construction and agronomic consistency.

[0025] In step S3, the judgment process for the preset similarity condition is as follows: calculate the cosine similarity of the soil attribute vectors of any two plots. When the cosine similarity is greater than the preset similarity threshold of 0.8, establish an attribute similarity edge between the two plot nodes. The soil attribute vector includes soil organic matter content, pH value, and permeability coefficient. By using soil organic matter, pH value, and permeability to form attribute vectors and establishing attribute similarity edges according to the cosine similarity threshold of 0.8, the geographical distance limitation can be overcome, effectively associating plots with highly consistent soil conditions. This fully explores the commonalities in soil productivity and similarities in growth patterns between non-adjacent plots, enabling the graph network to not only utilize spatial proximity information but also integrate soil attribute similarity, further improving the comprehensiveness and rationality of feature aggregation, effectively enhancing the feature complementarity capability in scenarios with missing data, and making the prediction model more consistent with the objective differences and inherent connections of farmland soil environment.

[0026] In step S4, a two-layer graph convolutional network is used for neighborhood feature aggregation. The specific calculation process is as follows: In the formula, The adjacency matrix A, containing physically adjacent edges and edges with similar attributes, is created by adding self-loops to the heterogeneous graph of the land parcels. Here, D is the degree matrix corresponding to the adjacency matrix A. The input features are the graph convolutional features of the l-th layer, where the input features are the graph convolutional features of the 0th layer. This represents the dynamic growth feature vector of the corresponding land parcel. Let be the trainable weight matrix of the l-th layer graph convolution. It is a ReLU nonlinear activation function. For the first The output features of the second-layer graph convolution are used as the neighborhood fusion growth features of the corresponding plots. The output feature dimension is 16. The two-layer graph convolution network performs neighborhood feature aggregation based on a heterogeneous graph containing physical adjacency edges and attribute similarity edges. It can weightedly fuse the dynamic growth features of the target plot itself with the effective features of adjacent and soil-similar plots, making full use of the spatial and environmental associations between plots to achieve information complementarity and feature enhancement. It effectively suppresses the feature loss problem caused by single plot data noise, cloud and fog occlusion, etc. At the same time, the two-layer graph convolution can fully extract the deep association features between plots while controlling the computational complexity, providing high-quality neighborhood fusion growth features with both spatial consistency and environmental similarity for subsequent yield prediction, significantly improving the robustness and prediction accuracy of the model.

[0027] In step S5, the calculation process of the basic productivity coefficient is as follows: based on each environmental factor in the static environmental data of the land parcel, each environmental factor is mapped to the [0,1] interval to obtain the fitness function value corresponding to each environmental factor; the minimum value among all fitness function values ​​is taken as the basic productivity coefficient of the corresponding land parcel. ,Right now: In the formula, This represents the fitness function value corresponding to the soil organic matter content. For soil The basic productivity coefficient is calculated by taking the minimum fit function value for each environmental factor. Strictly adhering to the principle of the agronomic short-board effect, it can determine the overall basic productivity of a plot based on the most limited soil environmental factors, avoiding the overestimation of the crop's potential due to short-term good growth in purely data-driven models. It embeds static environmental knowledge such as soil organic matter and pH value into the model in an explicit physical constraint manner, making the yield prediction results more consistent with the actual soil fertility of farmland, improving the scientificity, rationality and interpretability of the prediction, correcting characteristic biases at the mechanistic level, and ensuring that the final yield results do not violate agronomic common sense.

[0028] The fitness function corresponding to the soil organic matter content is: In the formula, The soil organic matter content of the plot. , This is a constant set based on agronomic experience in the target region, with a typical value range of [missing value]. , The soil The fitness function corresponding to the value is: In the formula, Soil of the plot value, For optimal crop growth value, This is a function sensitivity control parameter, with a typical value range of [value range missing]. By using these two fitness functions to standardize the mapping of soil organic matter and pH value respectively, environmental factors with different dimensions and value ranges can be uniformly converted to the [0,1] interval, ensuring the comparability of each factor in productivity evaluation. The organic matter adopts an S-shaped function to fit its gradual promoting law on crop growth, while the pH value adopts a Gaussian function to accurately match the crop's optimal growth response near neutral and its rapid decline after deviation. The parameters are set in combination with local agronomic experience, making the fitness calculation of environmental factors more consistent with the actual farmland mechanism. This provides a scientific, stable, and interpretable quantitative basis for subsequent short-board effect calculation and feature correction, significantly improving the agronomic rationality and physical consistency of the prediction results.

[0029] In step S5, the specific calculation method for constraining and correcting the neighborhood fusion growth characteristics of the corresponding plot using the basic productivity coefficient is as follows: In the formula, The final predicted features of the land parcel, To reflect the neighborhood integration and growth characteristics of the corresponding plots, This represents the basic productivity coefficient of the corresponding land parcel. To perform element-wise tensor product operations, the neighborhood fusion growth characteristics and basic productivity coefficients are constrained and corrected through element-wise tensor product. This allows for rigid calibration and dynamic suppression of remotely sensed growth characteristics based on the actual soil fertility upper limit of the plot. This avoids the overestimation of growth due to poor short-term growth and soil fertility caused by purely data-driven models. It organically combines the agronomic bottleneck mechanism with the spatial characteristics of graph networks, so that the final predicted features simultaneously take into account the dynamics of plot growth, spatial correlation information, and static soil productivity constraints. This preserves the discriminative power of growth characteristics and ensures that the prediction results conform to agronomic laws, significantly improving the reliability of features and the physical interpretability of model output.

[0030] The output regression model employs a random forest regression model, the configuration of which includes: 300-500 decision trees, a maximum depth of 15-20 layers, and a node splitting feature count of [missing information]. The split evaluation criterion uses mean squared error, and Bootstrap resampling is enabled for sample sampling. The model outputs the mean of all decision tree predictions as the crop yield per unit area for the corresponding plot. A random forest regression model with these parameters is used for yield prediction. This model leverages the ensemble learning advantage of 300-500 decision trees to reduce overfitting risk. A maximum depth of 15-20 layers controls complexity while ensuring model fitting ability. The node splitting feature number, mean square error splitting evaluation, and Bootstrap resampling mechanism can fully explore the nonlinear mapping relationship between growth, spatial correlation, and soil productivity in the final prediction features. The output multi-tree mean result can further smooth the prediction fluctuation of a single tree, improve the stability, accuracy, and generalization ability of plot-level yield estimation, adapt to the distribution characteristics of multi-source heterogeneous agricultural data, and meet the actual needs of refined management for prediction reliability.

[0031] This embodiment uses a region in my country as the area to be predicted, with winter wheat as the target crop. The prediction target is the yield per unit area (unit: kg / mu) of all winter wheat-grown plots within the county, to support the local agricultural and rural affairs bureau's work on refined winter wheat planting management, precise agricultural input allocation, grain storage planning, and disaster loss assessment. This embodiment fully discloses the entire process implementation details of the method. Those skilled in the art can repeat the technical solution of this invention without creative effort based on this embodiment and obtain the corresponding yield prediction results.

[0032] The basic data source and preprocessing method required for this embodiment are as follows: Multi-temporal remote sensing imagery was obtained from Sentinel-2 L2A atmospheric-corrected images with a spatial resolution of 10 meters. The temporal coverage spanned the entire growth cycle of winter wheat (from sowing in October of the previous year to maturity in June of the current year). Two cloudless or cloud-free images (less than 10%) were selected each month, resulting in a total of 16 valid temporal images. Image preprocessing included radiometric calibration, atmospheric correction, geometric fine correction (correction error less than 1 pixel), and study area cropping. All images were processed using both the WGS84 geographic coordinate system and the UTM projection coordinate system.

[0033] The land parcel boundary data uses vector boundary data of farmland parcels confirmed by the local agricultural and rural affairs department. The format is Shapefile polygon vector data, with each parcel corresponding to a unique parcel ID. The attribute table contains information such as the administrative village to which the parcel belongs, the confirmed area, and the planting type. The vector data is subjected to topological checks to eliminate topological errors such as overlap, gaps, and self-intersections, ensuring consistency with the spatial registration accuracy of the remote sensing image.

[0034] The static environmental data of the plots were obtained from the second national soil census and county-level farmland quality grade survey, with a spatial resolution of 30 meters. After spatial registration and plot-level statistics, soil attribute data for each plot were obtained, including soil organic matter content (unit: g / kg), soil pH value, and soil permeability coefficient (unit: cm / h). At the same time, the elevation and slope data of the plots were also obtained as auxiliary environmental factors.

[0035] The ground-based measured yield sample data uses the actual yield data of winter wheat plots harvested and measured by the local agricultural technology extension center for three consecutive years, accumulating 2,860 valid plot samples. Each sample includes information such as plot ID, measured yield per unit area, planting variety for the corresponding year, and water and fertilizer management level. Among them, 80% of the samples are used for model training, and 20% of the samples are used for model accuracy verification, ensuring that the samples cover plots with different soil fertility levels (high, medium, and low) and are spatially evenly distributed.

[0036] Detailed Implementation Steps of the Method Step S1: Extraction of time-series data of plot-level remote sensing statistical features and vegetation indices The core of this step is to transform pixel-level remote sensing data into plot-level features that conform to agricultural management logic, while eliminating interference from mixed pixels at plot edges. The specific implementation process is as follows: The core area of ​​the plot is screened by traversing the polygon of each confirmed plot. Using the GIS buffer analysis tool, a 5-meter negative buffer is generated inward for each plot polygon. Only the central area inside the plot that is more than 5 meters away from the boundary is retained. This removes mixed pixels that are at the edge of the plot and the ridges of the field or adjacent plots, ensuring that the extracted features only reflect the actual growth status of crops in the core planting area of ​​the plot.

[0037] Remote sensing statistical feature extraction targets the central area of ​​each plot, and statistically analyzes the statistical features of all effective pixels within the central area in the four core bands of the Sentinel-2 image on a time-phase basis. The four bands are: B2 band (blue light, B), B3 band (green light, G), B4 band (red light, R), and B8 band (near-infrared, NIR). The statistical indicators include the mean and standard deviation of the pixel values ​​in each band, which characterize the overall spectral reflectance level and growth uniformity of the plot.

[0038] Vegetation index time-series data calculation and extraction: For each scene of temporal image, the Normalized Difference Vegetation Index (NDVI) and Enhanced Vegetation Index (EVI) are calculated pixel by pixel. The calculation formula is as follows: , In the formula, NIR represents the near-infrared band pixel value, R represents the red band pixel value, and B represents the blue band pixel value. For the central area of ​​each plot, the average pixel values ​​of NDVI and EVI are calculated for each time phase to form the NDVI time series data and EVI time series data corresponding to each plot. The time series length is consistent with the number of time phases in the remote sensing image (16 time phases in this embodiment).

[0039] Texture feature extraction is performed on the NDVI data of each temporal image. Gray-Level Co-occurrence Matrix (GLCM) is used to calculate texture features with a 3×3 pixel calculation window and a step size of 1. This parameter combination preserves texture details within the plot while suppressing single-pixel noise interference. Two texture indices, entropy and contrast, are calculated for each pixel. For the central region of each plot, the pixel mean of the temporal texture indices is statistically analyzed and incorporated into the plot's remote sensing statistical feature system. Finally, all the above statistical features are arranged in temporal order to form the time-series data of remote sensing statistical features and vegetation indices corresponding to each plot, serving as the basic input for subsequent steps.

[0040] Step S2: Reconstruction of phenological time-series features and generation of dynamic growth feature vectors The core of this step is to extract agronomically significant crop phenological characteristics from NDVI time-series data, and to accurately quantify the growth dynamics of winter wheat throughout its entire growth cycle. The specific implementation process is as follows: For the NDVI time series data of each plot, Savitzky-Golay filtering (SG filtering) was used for smoothing. The filter window was set to 5 and the polynomial order to 2. This configuration can effectively smooth noise while preserving the phenological inflection point characteristics of the NDVI time series curve. The filter window was set to 5 and the polynomial order was set to 2 to eliminate time series fluctuations caused by cloud cover, sudden changes in light intensity, and sensor noise. At the same time, the 3σ criterion was used to remove abnormally low values ​​in the time series. That is, when the NDVI value of a certain time phase is lower than the time series mean minus 3 times the standard deviation, the NDVI mean of the two adjacent time phases was used for interpolation and replacement to ensure that the time series curve fits the actual growth and development pattern of winter wheat.

[0041] Phenological characteristic index extraction is based on the smoothed NDVI time series curve, extracting core phenological characteristic indicators for winter wheat, specifically including: The time of regrowth: the Julian day corresponding to the inflection point of the NDVI time series curve that marks a continuous rise from the trough of the overwintering period; The time of the jointing stage: the Julian day corresponding to the peak rate of the NDVI curve's rise; Peak at heading stage: the maximum value of the NDVI time-series curve; Heading period: Julian day corresponding to the peak of NDVI; NDVI decay rate during grain filling stage: the average daily rate of decrease in NDVI value from heading stage to maturity; Cumulative NDVI over the entire reproductive period: the area under the smoothed NDVI time-series curve.

[0042] The dynamic growth feature vector generation process involves min-max standardization of the six phenological period characteristic indicators, mapping them to the [0,1] interval, and arranging them in a fixed order to form a 6-dimensional dynamic growth feature vector for each plot, which serves as the node input feature for subsequent graph convolutional networks. The standardization formula is as follows: In the formula, For the original value of the feature, This is the minimum value of this feature among all land parcel samples. This is the maximum value of this feature across all land parcel samples.

[0043] Step S3: Construction of Heterogeneous Plot Graph Driven by Dual-Type Edges This step uses each parcel as a graph node to construct a heterogeneous parcel graph that includes both physical adjacency edges and attribute similarity edges, breaking down information silos between individual parcels. The specific implementation process is as follows: The graph node definition defines each land parcel in the area to be predicted as a node in the heterogeneous graph. The total number of nodes is equal to the total number of land parcels in the area to be predicted. Each node corresponds to a unique land parcel ID. The initial features of the nodes are the dynamically grown feature vectors generated in step S2.

[0044] Physical adjacency edges are constructed to connect land parcel nodes whose geographical locations meet preset adjacency conditions. This embodiment sets two-dimensional preset adjacency conditions, and a physical adjacency edge can be established if either condition is met: Condition 1: The two plots are directly adjacent in vector space, that is, the polygons of the two plots have a common edge; Condition 2: The centroid geographical Euclidean distance between two plots is less than a preset distance threshold of 100 meters. For all plot node pairs that meet the condition, an undirected physical adjacency edge is established between the nodes, with the edge weight uniformly set to 1.

[0045] Attribute similarity edges are constructed to connect land parcel nodes whose soil attributes meet preset similarity conditions. The specific implementation process is as follows: Construct a soil attribute vector for each plot, with a dimension of 3, arranged in the order of [soil organic matter content, soil pH value, soil permeability coefficient], and perform min-max standardization on each dimension to map it to the [0,1] interval; Traverse all pairwise plot nodes and calculate the cosine similarity of the standardized soil attribute vectors between the two plots. The calculation formula is as follows: In the formula, Let be the value of the i-th dimension of the soil attribute vector for plot A. Let be the value of the i-th dimension of the soil attribute vector for plot B, and n be the vector dimension (n=3 in this embodiment). The preset similarity threshold is set to 0.8. When the cosine similarity between two plots is greater than 0.8, an undirected attribute similarity edge is established between the two plot nodes, and the edge weight is uniformly set to 1.

[0046] The adjacency matrix generation of the heterogeneous graph is based on the aforementioned physical adjacency edges and attribute similarity edges. An adjacency matrix A is generated for the heterogeneous graph of the land parcels, with a dimension of N×N, where N is the total number of land parcel nodes. If two nodes have a physical adjacency edge or an attribute similarity edge, the corresponding element in the matrix has a value of 1; otherwise, it has a value of 0. Simultaneously, an identity matrix I is added to the adjacency matrix to introduce node self-loops, preventing the loss of node features during graph convolution. The final adjacency matrix is ​​then obtained for computation. .

[0047] Step S4: Neighborhood feature aggregation processing based on graph convolutional networks This step employs a two-layer graph convolutional network (GCN) to perform neighborhood aggregation on node features based on the heterogeneous parcel graph, thereby enhancing and completing the parcel features. The specific implementation process is as follows: The core calculation rules for graph convolutional networks: The forward propagation calculation formula for the graph convolutional layer used in this embodiment is as follows: The complete explanation of each parameter in the formula is as follows: The adjacency matrix of the heterogeneous plot graph after adding self-loops has a dimension of N×N; Adjacency matrix The corresponding degree matrix is ​​a diagonal matrix with dimensions N×N, where the value of the i-th element on the diagonal is... The sum of all elements in the i-th row; The input feature matrix of the l-th layer graph convolution has dimensions of Where N is the total number of nodes. The feature dimension of layer l; the input feature of layer 0. The matrix formed by the dynamic growth feature vectors of all plots, with dimension . ; The trainable weight matrix of the l-th layer graph convolution, with dimension . ,in Let l be the output feature dimension of the l-th layer; : Non-linear activation function; in this embodiment, the ReLU activation function is used, and its expression is: ; : No. The output feature matrix of the layer graph convolution has a dimension of .

[0048] This embodiment uses a two-layer graph convolutional network to perform feature aggregation, with the specific configuration as follows: First layer graph convolution: Input feature dimension Output feature dimension Weight matrix The dimension is 6×32, the weights are initialized using a uniform Xavier distribution, and the activation function is ReLU; Second layer graph convolution: Input feature dimension Output feature dimension Weight matrix The dimension is 32×16, and the weights are initialized using a uniform Xavier distribution without setting an activation function. Finally, the output feature matrix of the second layer graph convolution is... As the neighborhood fusion growth characteristics of each plot, each plot corresponds to a 16-dimensional neighborhood fusion growth feature vector.

[0049] The graph convolutional network was pre-trained using land parcel data corresponding to the measured yield samples. The training process was conducted end-to-end and jointly optimized with the subsequent yield regression model. The Adam optimizer was used, with a learning rate of 0.001, a batch size of 64, and 200 training epochs. The mean squared error (MSE) loss function was used. Training was terminated early when the validation set loss no longer decreased for 20 consecutive epochs to avoid overfitting.

[0050] Step S5: Feature constraint correction based on the agronomic short-board effect This step calculates the basic productivity coefficient based on the static environmental data of the plot, and applies explicit agronomic mechanism constraints to correct the neighborhood fusion growth characteristics, so as to avoid the prediction results from contradicting common sense in agricultural production. The specific implementation process is as follows: The environmental factor fitness function is constructed and calculated for the core environmental factors that determine the basic productivity of a plot. Fitness functions are constructed for each factor, and the original values ​​of each factor are mapped to the interval [0,1]. The closer the fitness value is to 1, the better the fit of the factor to crop growth. The details are as follows: The soil organic matter content fitness function adopts an S-shaped growth curve function, which conforms to the gradual promoting effect of soil organic matter on winter wheat growth. The function expression is as follows: In the formula: The values ​​represent the measured soil organic matter content of the corresponding plots, in g / kg. The slope of the function is controlled by a parameter set to 0.8 in this embodiment for winter wheat in the Huang-Huai-Hai region. The inflection point value of the function represents the critical value of organic matter content suitable for winter wheat growth; in this example, it is set to 15 g / kg. The output value range of this function is... The higher the organic matter content, the closer the compatibility value is to 1.

[0051] The soil pH fit function adopts a Gaussian distribution function to accurately match the pattern of optimal growth of winter wheat in neutral soil environment and the rapid decline in growth fit after deviating from the optimal value. The function expression is as follows: In the formula: The values ​​are the measured soil pH values ​​for the corresponding plots. To determine the optimal pH value for winter wheat growth, this example sets it to 6.5. The sensitivity control parameter for the function is set to 1.2 in this embodiment; the output value range of the function is... The closer the pH value is to 6.5, the closer the compatibility value is to 1.

[0052] The soil permeability coefficient fit function is a piecewise linear function, adapted to the suitable range of soil permeability for winter wheat. The function expression is as follows: In the formula, The measured soil permeability coefficient for the corresponding plot is given in cm / h. The function output range is [0,1], and the fit reaches its maximum value of 1 when the permeability coefficient is greater than 3.5 cm / h. The calculation of the basic productivity coefficient strictly follows the agronomic short-plate effect principle, that is, the upper limit of the basic productivity of the plot is determined by the most limiting environmental factor. The minimum value among the fit function values ​​of all environmental factors is taken as the basic productivity coefficient of the corresponding plot. The calculation formula is as follows: In the formula, This represents the soil organic matter content suitability value. This refers to the soil pH suitability value. This represents the soil permeability coefficient fit value. The value range is [0,1]. The closer the value is to 1, the higher the upper limit of the basic productivity of the land parcel.

[0053] The feature constraint correction calculation applies element-wise constraints to the neighborhood fusion growth characteristics of the corresponding plot using the basic productivity coefficient, resulting in the final predicted characteristics of the plot. The specific calculation formula is as follows: In the formula: This is the final predicted feature vector of the land parcel, with dimensions consistent with the neighborhood fusion growth feature (16 dimensions in this embodiment). This is the neighborhood fusion growth feature vector of the corresponding plot, which is the output feature of the second layer graph convolution; This represents the basic productivity coefficient of the corresponding land parcel; This involves element-wise tensor product operations, where each element of the neighborhood fusion growth feature vector is multiplied by the basic productivity coefficient to complete the global constraint correction of the feature.

[0054] Step S6: Training the yield regression model and outputting plot-level yield predictions This step uses a random forest regression model to construct a nonlinear mapping relationship between the final predicted features and the yield per unit area of ​​the land parcel, completing model training and yield prediction. The specific implementation process is as follows: The model training dataset is constructed for sample plots with actual ground yield labels. The final prediction features generated in step S5 are extracted as model input features, and the actual yield per unit area of ​​the corresponding plots are used as model labels to construct training and validation datasets. The training set accounts for 80% of the total samples (2288 samples), and the validation set accounts for 20% of the total samples (572 samples), ensuring that there are no significant differences in soil fertility distribution and spatial distribution between the training and validation sets.

[0055] The specific configuration parameters of the random forest regression model used in this embodiment are as follows: Number of decision trees (n_estimators): Set to 400; Maximum tree depth (max_depth): Set to 18 levels; Minimum number of samples for node split (min_samples_split): Set to 2; Minimum number of samples per leaf node (min_samples_leaf): Set to 1; Number of node split features (max_features): Set to log2M, where M is the dimension of the input features (in this example, M=16, corresponding to max_features=4). Split evaluation criteria: Mean squared error (MSE) is used; Sample sampling: Enable Bootstrap resampling and enable out-of-bag (OOB) scoring; Parallel computing: Set to 8 threads to improve training efficiency.

[0056] Model training and accuracy validation: The random forest regression model is trained using the training dataset, and its accuracy is validated using the validation dataset after training. The core validation metrics include: Coefficient of determination : Measures the model's ability to explain yield variation; the closer to 1, the better the fit. Root Mean Square Error (RMSE): Measures the overall deviation between predicted and measured values, expressed in kg / acre. Mean Absolute Percentage Error (MAPE): Measures the relative error of the prediction, expressed as a percentage. The model trained in this embodiment, and the validation set... The accuracy reached 0.92, the RMSE was 32.6 kg / mu, and the MAPE was 7.2%, meeting the precision requirements of agricultural precision management.

[0057] The plot-level yield prediction output targets all plots within the prediction area. The final prediction features generated in step S5 are input into a pre-trained random forest regression model. The model outputs the mean of all decision tree prediction results as the winter wheat yield per unit area prediction result for the corresponding plot, in kg / mu. Finally, the prediction results are associated with the plot ID and plot vector boundary to generate plot vector result data with yield attributes. A plot-level yield prediction statistical table is output simultaneously to support various decisions in precision agricultural management.

Claims

1. An image recognition and crop yield prediction method for precision agricultural management, characterized in that, Includes the following steps: S1 uses multi-temporal remote sensing images of the area to be predicted, combined with pre-acquired plot boundary data, to extract remote sensing statistical features and vegetation index time series data corresponding to each plot. S2 reconstructs the phenological time series characteristics of each plot based on the vegetation index time series data, and generates the dynamic growth feature vector of the corresponding plot. S3 uses each plot as a graph node to construct a plot heterogeneous graph; the connecting edges of the plot heterogeneous graph include physical adjacency edges and attribute similarity edges, wherein the physical adjacency edges are used to connect plot nodes whose geographical locations meet preset adjacency conditions, and the attribute similarity edges are used to connect plot nodes whose soil attributes meet preset similarity conditions. S4 uses a graph convolutional network to perform neighborhood feature aggregation on the dynamic growth feature vectors of each plot node based on the plot heterogeneous map, and obtains the neighborhood fusion growth features of each plot. S5. Based on the static environmental data of the plots in the area to be predicted, calculate the basic productivity coefficient of each plot, and use the basic productivity coefficient to constrain and correct the neighborhood fusion growth characteristics of the corresponding plots to obtain the final prediction characteristics of each plot; S6. Input the final prediction characteristics into the pre-trained yield regression model and output the crop yield per unit area prediction results for each plot.

2. The image recognition and crop yield prediction method for precision agricultural management according to claim 1, characterized in that, In step S1, the specific process of extracting the remote sensing statistical features and vegetation index time series data corresponding to each plot is as follows: traverse each plot polygon, select only the central area pixels inside the plot that are more than 5 meters away from the plot boundary, and statistically analyze the spectral mean and standard deviation of the central area pixels in the R, G, B, and NIR bands of the multi-temporal remote sensing images, as well as the mean values ​​of NDVI and EVI vegetation indices, and the entropy and contrast texture features calculated based on the gray-level co-occurrence matrix. The statistical average value is used as the remote sensing statistical features and vegetation index time series data of the corresponding plot.

3. The image recognition and crop yield prediction method for precision agricultural management according to claim 1, characterized in that, The specific process of S2 is as follows: the NDVI vegetation index time series data of each plot is filtered to remove abnormally low values; based on the smoothed NDVI time series curve, phenological period characteristic indicators including the greening period time, jointing period time, heading period peak, heading period time, grain filling period NDVI decay rate, and cumulative NDVI value throughout the entire growth period are extracted, and the phenological period characteristic indicators are used to construct the dynamic growth characteristic vector of the corresponding plot.

4. The image recognition and crop yield prediction method for precision agricultural management according to claim 1, characterized in that, In S3, the preset adjacency condition is: the two plots are directly adjacent in geographical location, or the centroid geographical Euclidean distance between the two plots is less than a preset distance threshold of 50-200 meters.

5. The image recognition and crop yield prediction method for precision agricultural management according to claim 1, characterized in that, In S3, the judgment process of the preset similarity condition is as follows: calculate the cosine similarity of the soil attribute vectors of any two plots, and when the cosine similarity is greater than the preset similarity threshold of 0.8, establish an attribute similarity edge between the two plot nodes; wherein, the soil attribute vector includes soil organic matter content, pH value, and permeability coefficient.

6. The image recognition and crop yield prediction method for precision agricultural management according to claim 1, characterized in that, In step S4, a two-layer graph convolutional network is used for neighborhood feature aggregation. The specific calculation process is as follows: In the formula, Adding self-loops to the heterogeneous graph of the land parcels results in an adjacency matrix containing physically adjacent edges and edges with similar attributes, where D is the degree matrix corresponding to the adjacency matrix A. The input features are the graph convolutional features of the l-th layer, where the input features are the graph convolutional features of the 0th layer. This represents the dynamic growth feature vector of the corresponding land parcel. Let be the trainable weight matrix of the l-th layer graph convolution. It is a ReLU nonlinear activation function. For the first The output features of the second-layer graph convolution are used as the neighborhood fusion growth features of the corresponding plots. The output feature dimension is 16-32. The adjacency matrix A is the adjacency matrix after adding node self-loops, i.e. ,in This is the adjacency matrix composed of the original physical adjacent edges and edges with similar attributes. It is an identity matrix.

7. The image recognition and crop yield prediction method for precision agricultural management according to claim 1, characterized in that, In step S5, the calculation process of the basic productivity coefficient is as follows: based on each environmental factor in the static environmental data of the land parcel, each environmental factor is mapped to the [0,1] interval to obtain the fitness function value corresponding to each environmental factor; the minimum value among all fitness function values ​​is taken as the basic productivity coefficient of the corresponding land parcel. ,Right now: In the formula, This represents the fitness function value corresponding to the soil organic matter content. For soil The corresponding fitness function value. This represents the fit function value corresponding to the soil permeability coefficient.

8. The image recognition and crop yield prediction method for precision agricultural management according to claim 7, characterized in that, The fitness function corresponding to the soil organic matter content is: In the formula, The soil organic matter content of the plot. , This is a constant set based on agronomic experience in the target region, with a typical value range of [value range missing]. ; , The standard unit of measurement for soil organic matter content is expressed as the number of grams of organic matter per kilogram of soil; the soil... The fitness function corresponding to the value is: In the formula, Soil of the plot value, For optimal crop growth value, This is a function sensitivity control parameter, with a typical value range of [value range missing]. .

9. The image recognition and crop yield prediction method for precision agricultural management according to claim 1, characterized in that, In step S5, the specific calculation method for constraining and correcting the neighborhood fusion growth characteristics of the corresponding plot using the basic productivity coefficient is as follows: In the formula, The final predicted features of the land parcel, To reflect the neighborhood integration and growth characteristics of the corresponding plots, This represents the basic productivity coefficient of the corresponding land parcel. This is an element-wise tensor product operation.

10. The image recognition and crop yield prediction method for precision agricultural management according to claim 1, characterized in that, The output regression model employs a random forest regression model, configured as follows: the number of decision trees is set to 300-500, the maximum depth is set to 15-20 layers, and the number of node split features is set to the input feature dimension. The split evaluation criterion adopts the mean squared error, and Bootstrap resampling is enabled for sample sampling; the model outputs the mean of all decision tree prediction results as the crop yield per unit area for the corresponding plot.

Images