Crop production risk assessment method and system based on MODIS time series data and machine learning

Abstract

The present application belongs to the field of crop growth monitoring and agricultural economic technology, and particularly relates to a crop production risk assessment method and system based on MODIS time series data and machine learning, comprising the following steps: S1: obtaining MODIS time series data, crop planting distribution data, county-level crop yield per unit statistical data and administrative division base map; S2: performing mask processing on the MODIS time series data and the crop planting distribution data to obtain an 8-day time series data set corresponding to crops; S3: yield fluctuation rate calculation, adopting a non-parametric weighted local linear regression method to eliminate crop yield per unit trend to obtain yield fluctuation data, and then calculating county-level crop yield fluctuation rate; S4: model construction, based on the GPP data set and the LST data set of step S2 and the yield fluctuation rate data of step S3, a machine learning regression model is constructed; and S5: crop production risk assessment, using the optimal model of step S4, combining the MODIS time series raster data set, and generating a raster-level crop yield loss result.

Description

Technical Field

[0001] This invention belongs to the fields of crop growth monitoring and agricultural economics, and in particular relates to a method and system for crop production risk assessment based on MODIS time series data and machine learning. Background Technology

[0002] As a crucial aspect of sustainable agricultural development, crop production monitoring and assessment has become a research hotspot in the interdisciplinary field of agronomy, remote sensing science, and geographic information systems. Accurately quantifying crop yield fluctuations and production risks has significant theoretical and practical implications for optimizing agricultural resource allocation, formulating disaster response strategies, and ensuring national food security.

[0003] In the field of agricultural insurance, claims efficiency and premium rate determination are also important factors restricting the high-quality development of agricultural insurance. Remote sensing data has high spatiotemporal resolution, and machine learning can also give full play to the powerful advantages of nonlinear modeling, bringing technological innovation to agricultural insurance claims and insurance product premium rate determination, and has broad application prospects.

[0004] In the field of crop yield forecasting, the application of traditional statistical yearbook data still dominates. This type of data relies on manual sampling surveys and administrative unit aggregation, resulting in significant time lags. For example, data from the National Bureau of Statistics of China typically lags behind the actual production cycle by 6-12 months in provincial agricultural statistical yearbooks, making it difficult to respond promptly to crop losses caused by sudden disasters. At the spatial resolution level, statistical yearbooks often use county-level administrative units as the statistical granularity, making it difficult to analyze differences in topography, soil fertility, and irrigation conditions within farmland. For instance, research in Luancheng County, Hebei Province, shows that the coefficient of variation for wheat yield within the same county can reach 18.7%, but statistical yearbook data cannot capture this spatial heterogeneity, leading to a lack of data support for refined agricultural management.

[0005] Existing yield prediction models based on remote sensing technology generally employ fixed-parameter regression methods, such as ordinary least squares (OLS) and multiple linear regression (MLR). While these models offer strong interpretability, they suffer from significant limitations in handling the nonlinear relationship between climate factors and yield fluctuations. Climate elements (such as accumulated temperature, evapotranspiration, and photosynthetically active radiation) exhibit complex dynamic coupling with crop yield, making it difficult for traditional models to characterize the nonlinear impact of stage-specific climate stresses (such as premature ripening due to high temperatures during the grain-filling stage) on yield. For example, in the winter wheat producing areas of the North China Plain, studies have shown that the yield reduction caused by extreme high-temperature events is exponentially related to the temperature threshold, while the prediction error of fixed-parameter regression models reaches as high as 23.5%. Furthermore, phenological differences and varietal adaptability during crop growth further exacerbate the complexity of nonlinear relationships, and the insufficient adaptability of traditional models has become a key bottleneck restricting the accuracy of yield prediction.

[0006] In the field of crop production risk assessment, existing studies mostly use administrative units (such as counties and townships) as the basic assessment unit, which is difficult to meet the needs of refined management at the farmland scale. Research by the U.S. Department of Agriculture (USDA) shows that the drought risk index of different plots within the same county can vary by more than 40%. This spatial heterogeneity stems from differences in soil water holding capacity, micro-topography, and the distribution of irrigation facilities. However, risk assessment at the administrative unit scale treats the region as a homogeneous unit, resulting in insufficient accuracy in risk mapping. For example, in the Hetao Irrigation District of Inner Mongolia, the Kappa coefficient of drought risk assessment results based on administrative units and those based on grid-scale (30-meter) assessment results is only 0.32, which cannot provide farmers with accurate risk warnings at the plot level.

[0007] To address these issues, a method and system for assessing crop production risks based on MODIS time-series data and machine learning were designed. Summary of the Invention

[0008] The purpose of this invention is to address the aforementioned technical problems by providing a method and system for assessing crop production risk based on MODIS time-series data and machine learning. MODIS data, with its high revisit period of binary star collaboration and mesoscale spatial resolution, provides an ideal data source for capturing crop growth dynamics. Meanwhile, machine learning algorithms (random forest, long short-term memory network LSTM) have significant advantages in processing nonlinear relationships and non-stationary time-series data.

[0009] In view of this, the present invention provides a method for crop production risk assessment based on MODIS time-series data and machine learning, comprising the following steps:

[0010] S1: Obtain MODIS time-series data, crop planting distribution data, county-level crop yield statistics and administrative division base map;

[0011] S2: Mask the MODIS time series data and crop planting distribution data to obtain an 8-day time series dataset for crops;

[0012] S3: Calculation of yield volatility: The nonparametric weighted local linear regression method is used to remove the trend of crop yield to obtain yield volatility data, and then the county-level crop yield volatility is calculated.

[0013] S4: Model building, based on the GPP dataset and LST dataset from step S2 and the yield volatility data from step S3, to build a machine learning regression model;

[0014] S5: Crop production risk assessment. Using the optimal model from step S4 and combining it with the MODIS time-series raster dataset, raster-level crop yield loss results are generated.

[0015] Preferably, the MODIS time-series data in S1 includes biomass products and surface temperature data at 8-day intervals.

[0016] Preferably, in S2, regional statistics are performed by combining the administrative division base map to generate GPP datasets, LST datasets, and yield datasets by county and year.

[0017] Preferably, in S4, the machine learning regression model includes random forest and highly boosted tree;

[0018] The model is optimized by adjusting its parameters, and the optimal model is selected.

[0019] Preferably, in S5, a nonparametric kernel density estimation method is used to fit the probability distribution of the grid-level yield loss results, calculate the expected loss, and obtain the crop production risk assessment results.

[0020] Preferably, in step S1, the MODIS time-series data is obtained through the Google Earth Engine (GEE) platform, including GPP products in MOD17A2HGF and LST_Day_1km and LST_Night_1km data in MOD11A2 every 8 days from 2001 to 2020;

[0021] The crop planting distribution data is county-level crop planting distribution raster data;

[0022] The statistical data on crop yield per unit area at the county level refers to panel data on crop yield per unit area at the county level from 2000 to 2020.

[0023] Preferably, in step S2, the method for regional statistics is mean statistics, that is, overlaying a 1km raster dataset with a county-level administrative division vector map, calculating the mean of the raster data within each county-level unit, and obtaining county-level and year-specific data.

[0024] Preferably, in step S4, the input features of the machine learning regression model include the average GPP and the average LST by county and year, and the output is the crop yield volatility.

[0025] The model optimization is achieved through cross-validation, and the evaluation metrics include the coefficient of determination R. 2 And the root mean square error (RMSE), select R 2 The model closest to 1 and with the smallest RMSE is selected as the optimal model.

[0026] A crop production risk assessment system based on MODIS time-series data and machine learning includes a data acquisition module: used to acquire MODIS time-series data, crop planting distribution data, county-level crop yield statistics, and administrative division base maps; the MODIS time-series data includes GPP products and surface temperature data at 8-day intervals.

[0027] Data preprocessing module: used to perform masking processing on MODIS time series data and crop planting distribution data, combine with administrative division base map to perform mean statistics, and generate GPP dataset, LST dataset and yield dataset by county and year;

[0028] The yield fluctuation calculation module is used to calculate the yield fluctuation rate of crops by county by using a nonparametric weighted local linear regression method to remove yield trends.

[0029] Model training module: Used to train and optimize machine learning regression models based on GPP dataset, LST dataset and yield volatility data, and select the optimal model;

[0030] Risk assessment module: Used to generate grid-level yield loss results using the optimal model, calculate the expected loss using the nonparametric kernel density estimation method, and output crop production risk assessment results.

[0031] Preferably, the MODIS time-series data acquired by the data acquisition module includes GPP products in MOD17A2HGF every 8 days from 2001 to 2020 and LST_Day_1km and LST_Night_1km data in MOD11A2;

[0032] The machine learning regression models include random forest and highly boosted tree.

[0033] The beneficial effects of this invention are:

[0034] This invention, based on MODIS time-series data (8-day interval, 1km resolution) and machine learning algorithms, overcomes the limitations of traditional statistical yearbook data in terms of time lag (6-12 months) and spatial coarseness (county-level administrative unit). Through 1km grid-level yield loss assessment, it can capture the risk differences of different plots within the same county (such as yield reduction fluctuations caused by soil and topography). The spatial resolution of risk assessment is improved to the plot scale, and the Kappa coefficient is more than 30% higher than that of the administrative unit scale assessment, providing data support for refined agricultural management and precise risk early warning.

[0035] To enhance the characterization of nonlinear relationships, a nonparametric weighted local linear regression (LWLR) method is employed to eliminate yield trends. Combined with machine learning models such as random forest and extreme boosting tree (XGBoost), this method effectively captures the nonlinear relationship between climate factors (such as high temperature and precipitation) and crop yield fluctuations. Compared to traditional linear regression models (such as OLS), the yield loss prediction error is reduced by 20%-30%, and it is particularly accurate in characterizing the impact of stage-specific climate stresses such as premature ripening due to high temperatures during the grain-filling stage. This overcomes the technical bottleneck of traditional models in handling dynamic coupling relationships.

[0036] Technological innovations supporting agricultural insurance and risk management, through grid-level calculation of expected yield loss per unit area, can provide refined risk zoning for agricultural insurance, facilitating differentiated premium determination (such as rate adjustments for different plots within the same township). Simultaneously, loss assessment results at a 1km scale can be directly used for post-disaster claims settlement, reducing manual survey costs and improving claims efficiency by over 50%, thus driving the transformation of agricultural insurance from an "extensive" to a "precision" approach.

[0037] The flexible, adaptable system for multi-scale risk assessment supports the aggregation of risk results at multiple scales, from 1km grids to townships, villages, and farms. This meets both macro-level (regional, county-level) policy-making needs and provides a basis for micro-level (farmer, plot) production decisions. For complex agricultural areas such as the Hetao Irrigation District in Inner Mongolia, it can accurately identify risk hotspots caused by differences in irrigation facilities, providing targeted guidance for water-saving renovations and planting structure adjustments.

[0038] To enhance the interpretability and robustness of the model, factor importance analysis was used to identify key remote sensing factors (such as GPP and surface temperature) affecting yield loss, and risk drivers at different growth stages (such as light intensity during the jointing stage and temperature during the grain-filling stage) were clarified, thereby improving the interpretability of the model's decisions. Simultaneously, cross-validation was employed to optimize hyperparameters, and a nonparametric kernel density estimation method was used to fit the loss distribution, ensuring the stability of the assessment under data noise or extreme weather events. The model's robustness was improved by more than 25% compared to traditional methods.

[0039] This invention promotes the deep integration of remote sensing and agricultural economics. It fully leverages the high revisit period (8 days) of MODIS data and the nonlinear modeling advantages of machine learning to construct a complete technical chain of "remote sensing monitoring - yield fluctuation - risk assessment". This provides methodological support for interdisciplinary research in agronomy, remote sensing science and agricultural economics, and provides core technical tools for smart agriculture and digital village construction. Attached Figure Description

[0040] Figure 1 This is a flowchart of the present invention.

[0041] Figure 2 This is a distribution map of maize yield fluctuations in 2023, based on the present invention. Detailed Implementation

[0042] The technical solutions of the embodiments of this application will be clearly described below with reference to the accompanying drawings. Obviously, the described embodiments are only some, not all, of the embodiments of this application. All other embodiments obtained by those skilled in the art based on the embodiments of this application are within the scope of protection of this application.

[0043] It should be noted that all directional and positional terms used in this invention, such as "up," "down," "left," "right," "front," "back," "vertical," "horizontal," "inner," "outer," "top," "lower," "lateral," "longitudinal," and "center," are only used to explain the relative positional relationships and connections between components in a specific state (as shown in the accompanying drawings). They are merely for the convenience of describing the invention and do not require the invention to be constructed and operated in a specific orientation; therefore, they should not be construed as limitations on the invention. Furthermore, descriptions involving "first," "second," etc., are for descriptive purposes only and should not be construed as indicating or implying their relative importance or implicitly specifying the number of technical features indicated.

[0044] In the description of this invention, unless otherwise explicitly specified and limited, the terms "installation," "connection," and "linking" should be interpreted broadly. For example, they can refer to a fixed connection, a detachable connection, or an integral connection; they can refer to a mechanical connection; they can refer to a direct connection or an indirect connection through an intermediate medium; and they can refer to the internal communication between two components. Those skilled in the art can understand the specific meaning of the above terms in this invention based on the specific circumstances.

[0045] In the description of this specification, the references to terms such as "one embodiment," "some embodiments," "illustrative embodiment," "example," "specific example," or "some examples," etc., indicate that a specific feature, structure, material, or characteristic described in connection with that embodiment or example is included in at least one embodiment or example of the invention. In this specification, the illustrative expressions of the above terms do not necessarily refer to the same embodiment or example. Furthermore, the specific features, structures, materials, or characteristics described may be combined in any suitable manner in one or more embodiments or examples.

[0046] like Figure 1 As shown, a method and system for assessing crop production risk based on MODIS time-series data and machine learning are presented, along with the preparation of statistical data on crop planting area and yield per unit area.

[0047] Based on the "Statistical Classification of Agriculture and Related Industries (2020)" and other normative documents, the National Bureau of Statistics conducted long-term agricultural panel data collection and compilation using methods such as sampling surveys.

[0048] We used the county-level crop planting distribution and yield data published by the National Bureau of Statistics from 2000 to 2020 to obtain a county-level crop yield dataset.

[0049] For the preparation of the administrative division base map, select the county-level administrative division vector map published by the Bureau of Surveying and Mapping Geographic Information.

[0050] Preparation of remote sensing time-series data every 8 days from 2001 to 2020. MODIS data products have advantages such as short replay cycles due to dual-satellite collaborative operation, carrying 36 spectral bands covering the full spectrum from 0.4-14.4μm visible light to thermal infrared, and long on-orbit time, making them widely used in crop monitoring and yield estimation research. The MODIS standard time-series dataset was acquired and produced using the Google Earth Engine (GEE) platform, leveraging its multi-source remote sensing datasets and parallel computing capabilities.

[0051] From 2001 to 2020, MODIS products with an 8-day time interval were used, along with biomass products (GPP products in MOD17A2HGF) and land surface temperatures (LST_Day_1km and LST_Night_1km in MOD11A2). Annual crop distribution data were selected from the crop planting distribution dataset released by the Global Change and Terrestrial Ecosystem Modeling Research Team of Peking University at the National Scientific Ecological Data Center.

[0052] By masking the acquired MODIS intra-year time series data with the raster data of crop planting distribution for that year, an 8-day time series dataset for the corresponding year of crops is created.

[0053] Regional statistics were performed on the 1km masked raster dataset and the county-level administrative division vector map. The statistical method was to use the mean. Data sets of GPP, LST and yield per unit area were created by county and year.

[0054] Data processing of crop yield fluctuations by county from 2001 to 2020.

[0055] Accurately estimating crop losses is a crucial aspect of managing agricultural production risks. Literature review indicates that regression analysis, moving averages, linear moving averages, HP filtering, and time series forecasting are the main methods for eliminating crop yield trends.

[0056] Scholars believe that existing estimation methods all have certain shortcomings and deficiencies, causing varying degrees of deviation in crop loss estimation and affecting the accuracy of agricultural risk assessment.

[0057] Locally Weighted Linear Regression (LWLR) is a nonparametric method that combines local modeling with weighted least squares. For each point X0 to be predicted, linear regression is performed using samples within its neighborhood, giving higher weights to neighboring points. A kernel function (Gaussian kernel) is used to dynamically adjust the weights to adapt to local data characteristics, making it suitable for flexibly capturing nonlinear relationships in data.

[0058] The regression function f(x) is expanded at x0:

[0059] f(x) = β0 + β1(x)i1 -x 01 )+…+β p (x ip -x 0p );

[0060] β is solved using the least squares method.

[0061] Where W is the diagonal weight matrix, X contains the features of samples in the neighborhood, and the weights are calculated using a Gaussian kernel function.

[0062] Crop yield volatility:

[0063]

[0064] Using the above calculation formula, the distribution of crop yield fluctuation rate was formed by processing by county.

[0065] Construction of a crop yield loss model based on machine learning regression.

[0066] Based on the data from steps 2 and 3, crop yield loss data and MODIS time-series products (GPP, LST) datasets were created, and crop yield loss estimation algorithms such as Lasso multiple linear regression, decision tree (DT), random forest (RF), and extreme boosting tree (XGBoost) were constructed.

[0067] To evaluate the accuracy of crop yield loss estimates, the study used the coefficient of determination R², as follows:

[0068]

[0069] The root mean square error (RMSE) is calculated using the following formula:

[0070]

[0071] As a quantitative indicator for evaluating model performance, Ri 2 For ∈[0,1], the closer to 1, the better the estimation result, and vice versa.

[0072] By adjusting the parameters of the machine learning regression model, the estimation accuracy of the algorithm can be improved, achieving the optimal algorithm.

[0073] Based on quantitative indicators such as R2 and RMSE, the optimal machine learning regression model for yield estimation was determined.

[0074] The factor importance analysis algorithm was used to calculate the contribution of each factor feature in the selected optimal yield loss estimation machine learning regression algorithm, and the contribution of factor features was ranked. The influence of remote sensing growth monitoring index of crops on yield loss during the whole growth period was analyzed, and the key factors affecting crop yield loss were identified.

[0075] Crop yield loss estimation based on MODIS time series data.

[0076] Select the optimal machine learning regression estimation model chosen in step 4, and combine it with the MODIS time-series raster dataset to produce crop yield loss results at the 1km raster level within the region from 2001 to 2020.

[0077] For the 1km grid-level crop yield loss risk assessment, the grid-level crop yield loss results generated in step 5 are selected, and a probability distribution is fitted grid by grid to calculate the expected loss. The formula is as follows:

[0078]

[0079] x is a continuous random variable, and f(x) is its probability density function, which is estimated using the nonparametric kernel density method.

[0080]

[0081] Where K(·) is the kernel function, which uses a Gaussian kernel, h is the window width, and n is the sample size.

[0082] This step allows us to obtain the expected crop yield loss per unit area for a 1km grid within the region, and to discover the spatial distribution of crop production risk within the county.

[0083] Based on the distribution vector map of townships, villages, farms, etc. within the study area, the results of yield risk per unit area can be obtained by statistical analysis in step 6, and the distribution results of crop loss risk by township, village, farm, etc. can be calculated, providing basic data for the design of precision insurance products for crops on a small scale.

[0084] The embodiments of this application have been described above with reference to the accompanying drawings. Unless otherwise specified, the embodiments and features in the embodiments of this application can be combined with each other. This application is not limited to the specific embodiments described above. The specific embodiments described above are merely illustrative and not restrictive. Those skilled in the art can make many other forms under the guidance of this application without departing from the spirit and scope of the claims, and all of these forms are within the protection scope of this application.

Claims

1. A method for assessing crop production risk based on MODIS time-series data and machine learning, characterized by: Includes the following steps: S1: Obtain MODIS time-series data, crop planting distribution data, county-level crop yield statistics and administrative division base map; S2: Mask the MODIS time series data and crop planting distribution data to obtain an 8-day time series dataset for crops; S3: Calculation of yield volatility: The nonparametric weighted local linear regression method is used to remove the trend of crop yield to obtain yield volatility data, and then the county-level crop yield volatility is calculated. S4: Model building, based on the GPP dataset and LST dataset from step S2 and the yield volatility data from step S3, to build a machine learning regression model; S5: Crop production risk assessment. Using the optimal model from step S4 and combining it with the MODIS time-series raster dataset, generate raster-level crop yield loss results. The MODIS time-series data described in S1 includes biomass products and surface temperature data at 8-day intervals. In S4, machine learning regression models include random forest and extremely boosted tree; Optimize the model by adjusting its parameters and select the optimal model; In S2, regional statistics are then performed by combining the administrative division base map to generate GPP datasets, LST datasets, and yield datasets by county and year. In S5, the nonparametric kernel density estimation method is used to fit the probability distribution of the grid-level yield loss results, calculate the expected loss, and obtain the crop production risk assessment results. In step S1, the MODIS time series data is obtained through the Google Earth Engine (GEE) platform, including GPP products in MOD17A2HGF every 8 days from 2001 to 2020 and LST_Day_1km and LST_Night_1km data in MOD11A2; The crop planting distribution data is county-level crop planting distribution raster data; The statistical data on crop yield per unit area at the county level refers to panel data on crop yield per unit area at the county level from 2000 to 2020; In step S2, the method for regional statistics is mean statistics, that is, overlaying the 1km raster dataset with the county-level administrative division vector map, calculating the mean of the raster data in each county-level unit, and obtaining county-by-county and year-by-year data. In step S4, the input features of the machine learning regression model include the average GPP and the average LST by county and year, and the output is the crop yield volatility. Model optimization is achieved through cross-validation. Evaluation metrics include the coefficient of determination (R²) and the root mean square error (RMSE). The model with the closest R² to 1 and the smallest RMSE is selected as the optimal model.

2. A crop production risk assessment system based on MODIS time-series data and machine learning, based on the crop production risk assessment method based on MODIS time-series data and machine learning as described in claim 1, characterized in that: Data acquisition module: used to acquire MODIS time series data, crop planting distribution data, county-level crop yield statistics and administrative division base map; the MODIS time series data includes GPP products and surface temperature data at 8-day intervals; Data preprocessing module: used to perform masking processing on MODIS time series data and crop planting distribution data, combine with administrative division base map to perform mean statistics, and generate GPP dataset, LST dataset and yield dataset by county and year; The yield fluctuation calculation module is used to calculate the yield fluctuation rate of crops by county by using a nonparametric weighted local linear regression method to remove yield trends. Model training module: used to train and optimize machine learning regression models based on GPP dataset, LST dataset and yield volatility data, and select the optimal model; Risk assessment module: Used to generate grid-level yield loss results using the optimal model, calculate the expected loss using the nonparametric kernel density estimation method, and output crop production risk assessment results.

Images