A grid scale crop yield prediction method based on air quality mechanism

Abstract

This invention discloses a grid-scale crop yield prediction method based on air quality mechanisms, belonging to the field of agricultural remote sensing technology. The method includes the following steps: S1, dividing the study area into grids and collecting basic data, then preprocessing the basic data sequentially; S2, constructing a physiological regulation function based on the preprocessed basic data and determining stomatal conductance; S3, determining air quality mechanism indicators and pollutant mechanism indicators based on stomatal conductance; S4, determining the predicted crop yield based on the air quality mechanism indicators and pollutant mechanism indicators, and outputting the baseline yield, total loss, relative loss rate, pollutant contribution rate, and spatially aggregated total yield. The final prediction target variable of this invention is yield, suitable for large-scale yield estimation at the grid scale; it can be used for grid-scale crop yield estimation, air pollution-induced damage assessment, and risk analysis.

Description

Technical Field

[0001] This invention relates to the field of agricultural remote sensing technology, specifically to a grid-scale method for predicting crop yield based on the mechanism of air quality. Background Technology

[0002] Changes in air quality have a significant impact on crop growth, development, and yield. Different pollutants exhibit distinct pathways of action on crops. Ozone enters the plant through leaf stomata and induces oxidative stress; particulate matter alters the canopy radiation environment and leaf processes through extinction, scattering, and deposition; while sulfur dioxide and nitrogen dioxide affect leaf physiological state, soil environment, and nutrient supply through leaf absorption and deposition. Therefore, the impact of air pollution on crop yield is not solely due to changes in concentration, but is closely related to the crop's growth stage, leaf exchange capacity, and environmental regulation conditions.

[0003] Existing crop yield estimation technologies primarily use meteorological, remote sensing, soil, and management information as inputs, employing statistical or machine learning models to predict yield. While some methods have begun to incorporate air pollutant indicators, they typically only directly input pollutant concentrations as general features into the model, lacking quantification of the differential effects of pollutants and failing to reflect the transformation mechanism from external exposure to the actual intensity of crop impact. Furthermore, crop responses to air pollution exhibit distinct temporal characteristics. Crop sensitivities vary at different growth stages, and stomatal opening is modulated by factors such as radiation, temperature, air aridity, and soil moisture. Current technologies rarely incorporate process parameters such as phenological identification, physiological regulation functions, and stomatal conductance into pollution intensity calculations, resulting in a lack of effective linkage between air pollutant concentrations and actual crop damage, thus affecting the mechanistic support and interpretability of yield predictions. Moreover, in regional or national scale applications, air quality, meteorological, remote sensing, soil, and management data often have different spatiotemporal resolutions. Therefore, existing technologies generally suffer from insufficient mechanistic expression, inadequate spatiotemporal scale linkages, and weak capabilities in loss assessment and pollutant contribution analysis.

[0004] Therefore, it is necessary to propose a grid-scale crop yield prediction method based on the mechanism of air quality action. This method incorporates crop phenological identification, physiological regulation functions, stomatal conductance, and differentiated indicators of the effects of different pollutants into the crop yield estimation process. Furthermore, it combines multi-source data and machine learning models to achieve crop yield prediction, pollution damage assessment, and pollutant contribution analysis. Summary of the Invention

[0005] To address the problems in existing technologies such as insufficient quantification of the effects of different air pollutants, insufficient sensitivity expression of phenological stages, difficulty in unifying multi-source data at the raster scale, and insufficient support for yield prediction mechanisms, this invention proposes a raster-scale crop yield prediction method based on the mechanism of air quality action.

[0006] The technical solution of this invention is: a grid-scale crop yield prediction method based on the mechanism of air quality includes the following steps:

[0007] S1. Divide the study area into grids and collect basic data. Then, preprocess the basic data in sequence.

[0008] S2. Based on the preprocessed basic data, construct the physiological regulation function and determine the stomatal conductance;

[0009] S3. Determine air quality mechanism indicators and pollutant mechanism indicators based on porosity;

[0010] S4. Based on air quality mechanism indicators and pollutant mechanism indicators, determine the predicted crop yield and output the baseline yield, total loss, relative loss rate, pollutant contribution rate, and spatial aggregate total yield.

[0011] Furthermore, in S1, preprocessing includes temporal scale unification, spatial scale unification, and missing value imputation;

[0012] The expression for time-scale uniformity is:

[0013] ;

[0014] in, This represents the daily average value. These are the original variable values ​​on an hourly scale; The number of valid hours for the day;

[0015] The expression for spatial scale unification is:

[0016] ;

[0017] ;

[0018] in, These are the interpolated raster variable values; For the site Observed values; For grid With the site The weights; For grid With the site The distance; This is the distance attenuation coefficient; The total number of sites participating in the interpolation;

[0019] Missing value imputation is specifically as follows:

[0020] If there are missing dates or hours, the expression for filling them in is:

[0021] ;

[0022] in, The corrected variable values; The variable value is from the previous day; The variable value for the next day.

[0023] If consecutive missing values ​​exceed a preset threshold, the expression for filling them in is:

[0024] ;

[0025] in, This represents the number of neighboring grid cells; For grid The neighborhood set; For neighborhood grid The variable value.

[0026] Furthermore, S2 includes the following sub-steps:

[0027] S21. Based on the preprocessed basic data, determine several growth stages according to the target crop type and identify crop phenology;

[0028] S22. Based on crop phenology, construct physiological regulation functions, specifically including phenological regulation function, radiation regulation function, temperature regulation function, saturated water vapor pressure deficit regulation function, and soil moisture content regulation function;

[0029] S23. Determine stomatal conductance based on phenological regulation function, radiation regulation function, temperature regulation function, saturated water vapor pressure deficit regulation function, and soil moisture content regulation function.

[0030] Furthermore, in S21, the method for determining the growth stage is as follows: The cumulative effective accumulated temperature is calculated based on the sowing date and the crop base temperature; when the cumulative effective accumulated temperature reaches a threshold, the crop is determined to have entered the growth stage; the cumulative effective accumulated temperature... The expression is:

[0031] ;

[0032] ;

[0033] in, The average daily temperature; Base temperature for crops; The date sequence number for accumulated temperature; For grid ,years Next stage The starting date sequence number; For grid ,years Next stage End date sequence number; For the stage The corresponding set of dates;

[0034] In S22, the phenological regulation function The expression is:

[0035] ;

[0036] in, This refers to the serial number of the start date of the reproductive period; The date number corresponding to the peak physiological activity; This refers to the serial number of the end date of the reproductive period;

[0037] In S22, the radiation adjustment function The expression is:

[0038] ;

[0039] in, Photosynthetically active radiation; The light intensity is half-saturated.

[0040] In S22, the temperature adjustment function The expression is:

[0041] ;

[0042] in, The average daily temperature; The minimum suitable temperature for crop growth; The optimal temperature for crops; This is the optimal temperature for crop growth.

[0043] In S22, the saturated vapor pressure deficit adjustment function The expression is:

[0044] ;

[0045] in, This indicates a saturated vapor pressure deficit. The threshold for low dryness is [not specified]. The high dryness threshold;

[0046] In S22, the soil moisture content adjustment function The expression is:

[0047] ;

[0048] in, This refers to the soil moisture content in the root zone. This refers to the moisture content at the wilting point. Field holding capacity;

[0049] In S23, porosity The expression is:

[0050] ;

[0051] in, This represents the maximum porosity.

[0052] Furthermore, S3 includes the following sub-steps:

[0053] S31. Based on stomatal conductance, cumulative effect of ozone concentration exceeding the plant sensitivity threshold, weighted cumulative effect of high-concentration ozone, and effective ozone dose exceeding the flux threshold, the comprehensive ozone mechanism index is determined by weighting.

[0054] S32. Calculate the ozone photosynthesis inhibition index based on the comprehensive ozone mechanism index and the maximum stomatal conductance.

[0055] S33. Calculate the volume extinction coefficient of particulate matter based on the concentrations of PM2.5 and PM10.

[0056] S34. Calculate direct radiation based on the volume extinction coefficient of particulate matter;

[0057] S35. Calculate the scattered radiation;

[0058] S36. Calculate the particulate matter radiation regulation index based on direct radiation and scattered radiation;

[0059] S37. Calculate the particulate matter settling amount;

[0060] S38. Calculate the amount of sulfur dioxide absorbed by the leaves and the amount of sulfur deposition.

[0061] S39. Calculate the soil pH at the end of the stage based on the amount of sulfur deposition, and determine the soil acidification index.

[0062] S311. Calculate the leaf absorption dose of nitrogen dioxide and the amount of nitrogen deposition;

[0063] S312. The comprehensive ozone mechanism index, ozone photosynthesis inhibition index, particulate matter radiation regulation index, and particulate matter deposition amount are used as air quality mechanism indexes, and the leaf absorption dose of sulfur dioxide, sulfur deposition amount, soil acidification index, leaf absorption dose of nitrogen dioxide, and nitrogen deposition amount are used as pollutant mechanism indexes.

[0064] Furthermore, in S31, the cumulative effect of ozone concentration exceeding the plant sensitivity threshold... The expression is:

[0065] ;

[0066] in, This refers to the hourly ozone concentration. For grid ,years ,stage Daytime hourly gatherings;

[0067] In S31, the weighted cumulative effect of high concentrations of ozone The expression is:

[0068] ;

[0069] in, It is an exponential function;

[0070] In S31, the effective ozone dose after exceeding the flux threshold The expression is:

[0071] ;

[0072] in, Unit conversion factor; This refers to the ozone flux threshold. The step size is in hours; Pore ​​conductance;

[0073] In S31, the comprehensive ozone mechanism index The expression is:

[0074] ;

[0075] in, This is the standardized AOT40; This is the standardized W126; For standardization ; The weights of the standardized AOT40; W126 is the weighted and standardized version; For standardization The weights;

[0076] In S32, the ozone photosynthesis inhibition index The expression is:

[0077] ;

[0078] in, Maximum porosity;

[0079] In S33, the volume extinction coefficient of particulate matter The expression is:

[0080] ;

[0081] in, PM2.5 extinction efficiency; PM10 extinction efficiency; PM2.5 concentration; PM10 concentration;

[0082] In S34, direct radiation The expression is:

[0083] ;

[0084] in, Potential radiation; Equivalent atmospheric optical path;

[0085] In S35, scattered radiation The expression is:

[0086] ;

[0087] in, The particulate scattering coefficient;

[0088] In S36, the particulate matter radiation regulation index The expression is:

[0089] ;

[0090] in, The average scattered radiation of the stage; The average potential radiation for the stage; The average direct radiation during the phase;

[0091] In S37, particulate matter sedimentation The expression is:

[0092] ;

[0093] ;

[0094] in, The daily time step; This refers to particulate matter settling flux; PM2.5 deposition rate; PM10 deposition velocity; For the stage The corresponding set of dates;

[0095] In S38, the leaf absorption dose of sulfur dioxide The expression is:

[0096] ;

[0097] in, This refers to the concentration of sulfur dioxide.

[0098] In S38, sulfur deposition The expression is:

[0099] ;

[0100] in, This refers to the settling velocity of sulfur dioxide. This refers to the concentration of sulfur dioxide. For the stage The corresponding set of dates;

[0101] In S39, the soil pH at the end of the stage The expression is:

[0102] ;

[0103] in, This represents the soil pH at the end of the previous stage. The corresponding coefficient for acidification;

[0104] In S39, soil acidification index The expression is:

[0105] ;

[0106] in, For reference soil pH;

[0107] In S311, nitrogen deposition The expression is:

[0108] ;

[0109] in, This refers to the settling velocity of nitrogen dioxide. This refers to the concentration of nitrogen dioxide.

[0110] In S311, the leaf absorption dose of nitrogen dioxide The expression is:

[0111] ;

[0112] in, This represents the concentration of nitrogen dioxide.

[0113] Furthermore, S4 includes the following sub-steps:

[0114] S41. Construct feature vectors based on air quality mechanism indicators and pollutant mechanism indicators;

[0115] S42. Based on the feature vectors, construct a crop yield prediction model and determine the predicted crop yield;

[0116] S43. Supervised learning training of crop yield prediction models using historical samples;

[0117] S44. Set the air quality mechanism index to the pollution-free baseline level to obtain the baseline output;

[0118] S45. Calculate the total loss based on the crop yield prediction model trained under supervised learning and the baseline yield;

[0119] S46. Calculate the relative loss rate based on the total loss.

[0120] S47. Calculate the pollutant contribution rate;

[0121] S46. Calculate the total spatial aggregate yield based on the crop yield forecast;

[0122] S47. Output the predicted crop yield, baseline yield, total loss, relative loss rate, pollutant contribution rate, and spatial aggregate total yield as results.

[0123] Furthermore, in S41, the eigenvector The expression is:

[0124] ;

[0125] in, This represents the cumulative effect of ozone concentration exceeding the plant sensitivity threshold. This represents the weighted cumulative effect of high concentrations of ozone. The effective ozone dose after exceeding the flux threshold; As a comprehensive indicator of ozone mechanism; As an indicator of ozone inhibition of photosynthesis; As an indicator of particulate matter radiation regulation; This refers to particulate matter settling volume; This refers to the leaf absorption dose of sulfur dioxide; This refers to the amount of sulfur deposition. Indicators of soil acidification; This refers to the leaf absorption dose of sulfur dioxide; This refers to nitrogen deposition. The average temperature of the stage; This refers to the cumulative precipitation over a given period. The average photosynthetically active radiation (APRA) during the phase; This represents the average soil moisture content for that stage. The leaf area index is the stage index; For stage NDVI; For staged solar-induced chlorophyll fluorescence; This is a soil attribute vector; To manage feature vectors;

[0126] In S42, crop yield forecasting The expression is:

[0127] ;

[0128] in, For machine learning model functions;

[0129] In S43, the expression for the objective function for supervised learning training is:

[0130] ;

[0131] in, For training the objective function; For the set of model parameters; This represents the actual crop yield; The regularization coefficient is used. It is a regularization term;

[0132] In S44, the baseline production The expression is:

[0133] ;

[0134] in, These are the feature vectors under no-pollution or baseline-pollution scenarios;

[0135] In S45, the total loss The expression is:

[0136] ;

[0137] In S46, the relative loss rate The expression is:

[0138] ;

[0139] In S47, the contribution rate of pollutants The expression is:

[0140] ;

[0141] in, Index for pollutant categories; pollutants The absolute contribution value in production forecasting;

[0142] In S46, the total output of spatial aggregation The expression is:

[0143] ;

[0144] in, For target crops in the grid ,years The planting area.

[0145] The beneficial effects of this invention are:

[0146] (1) This invention constructs intermediate indicators based on the differentiated action mechanisms of ozone, particulate matter, sulfur dioxide and nitrogen dioxide, rather than directly using pollutants as ordinary explanatory variables;

[0147] (2) The present invention incorporates phenological stages, stomatal conductance, radiation process, sedimentation process, and soil process into the algorithm flow;

[0148] (3) This invention uses mechanistic formulas for the parts of the mechanism that are clearly defined, and machine learning modeling for the parts of the mechanism that are not fully defined;

[0149] (4) The final prediction target variable of this invention is yield, which is suitable for large-scale yield estimation at the grid scale; it can be used for crop yield estimation at the grid scale, air pollution damage assessment and risk analysis. Attached Figure Description

[0150] Figure 1 This is a flowchart of a grid-scale crop yield prediction method based on the mechanism of air quality. Detailed Implementation

[0151] The embodiments of the present invention will be further described below with reference to the accompanying drawings.

[0152] like Figure 1 As shown, this invention provides a grid-scale crop yield prediction method based on the mechanism of air quality, comprising the following steps:

[0153] S1. Divide the study area into grids and collect basic data. Then, preprocess the basic data in sequence.

[0154] S2. Based on the preprocessed basic data, construct the physiological regulation function and determine the stomatal conductance;

[0155] S3. Determine air quality mechanism indicators and pollutant mechanism indicators based on porosity;

[0156] S4. Based on air quality mechanism indicators and pollutant mechanism indicators, determine the predicted crop yield and output the baseline yield, total loss, relative loss rate, pollutant contribution rate, and spatial aggregate total yield.

[0157] Step S1 is used to construct a unified research spatial unit and temporal index, and to collect, quality control, spatiotemporal registration and missing data of air quality, meteorology, remote sensing, soil, management and yield, forming a rasterized basic database that can be used for subsequent calculations;

[0158] Step S2 is used to identify the growth stage of the target crop in each grid based on sowing information, accumulated temperature process and remote sensing growth information, and to calculate the phenological regulation function, radiation regulation function, temperature regulation function, saturated vapor pressure deficit regulation function and soil moisture content regulation function related to stomatal opening, leaf activity and environmental stress.

[0159] Step S3 is used to establish differentiated mechanism indicators for ozone, particulate matter, sulfur dioxide and nitrogen dioxide respectively. Among them, ozone focuses on high-concentration cumulative exposure and stomatal flux mechanism, particulate matter focuses on radiation modulation and deposition mechanism, and sulfur dioxide and nitrogen dioxide focus on leaf absorption, deposition and soil or nutrient effect mechanism. The hourly or daily scale mechanism quantities are aggregated into growth stage scale characteristics.

[0160] Step S4 is used to construct a yield prediction feature vector by combining mechanistic indicators with meteorological, remote sensing, soil and management information, and use machine learning models to obtain the predicted yield of the target crop at the grid scale; at the same time, a pollution-free or baseline pollution scenario is constructed, and the baseline yield, loss amount, loss rate, pollutant contribution rate and administrative unit aggregation results are calculated.

[0161] In this embodiment of the invention, S1 includes preprocessing including time scale unification, spatial scale unification, and missing value imputation;

[0162] This step is used to establish a unified raster research unit, time index, and basic database to provide consistent data input for subsequent phenological identification, mechanistic index calculation, and yield forecasting.

[0163] Based on the study area and target resolution, the area is divided into regular grid cells. Each grid cell serves as a unified unit for subsequent air quality, meteorological, remote sensing, soil, management, and yield data. The study area is divided into... One grid cell, Number the grid cells. .

[0164] Establish a unified time index for different data types. Let... Number the year. To represent the total number of years studied, ; For the first The day number in the middle of the year, For the first Total number of days in a year ; For hourly serial numbers, ; Numbering the stages of crop growth. This refers to the total number of crop growth stages. .

[0165] The following data was collected from each grid cell:

[0166] Air quality data: ozone concentration, PM2.5 concentration, PM10 concentration, sulfur dioxide concentration, nitrogen dioxide concentration;

[0167] Meteorological data: temperature, photosynthetically active radiation, total radiation, precipitation, saturated vapor pressure deficit, soil moisture content;

[0168] Remote sensing data: NDVI, LAI, SIF;

[0169] Soil data: soil pH, organic matter, total nitrogen, available phosphorus, and available potassium;

[0170] Management data: sowing period, irrigation, fertilization, variety;

[0171] Production data, historical statistics on yield per unit area and total output.

[0172] For hourly-scale raw data, the hourly series is retained for calculating ozone exposure, stomatal conductance, and deposition flux. For daily-scale model input variables, if the original values ​​of the variables are hourly-scale, the daily average is calculated.

[0173] The expression for time-scale uniformity is:

[0174] ;

[0175] in, This represents the daily average value. These are the original variable values ​​on an hourly scale; The number of valid hours for the day;

[0176] Map station observation data, vector data, or raster data of different spatial resolutions to the target raster.

[0177] The expression for spatial scale unification is:

[0178] ;

[0179] ;

[0180] in, These are the interpolated raster variable values; For the site Observed values; For grid With the site The weights; For grid With the site The distance; This is the distance attenuation coefficient; The total number of sites participating in the interpolation;

[0181] Missing value imputation is specifically as follows:

[0182] Quality control is performed on missing or abnormal data to form a basic database.

[0183] When a single day or hour is missing, it can be filled using the average of adjacent time periods, the average of adjacent dates, or time interpolation. When consecutive missing values ​​exceed a preset threshold, it can be filled by combining neighboring grids, neighboring stations, or the average of the same period in history.

[0184] If there are missing dates or hours, the expression for filling them in is:

[0185] ;

[0186] in, The corrected variable values; The variable value is from the previous day; The variable value for the next day.

[0187] If consecutive missing values ​​exceed a preset threshold, the expression for filling them in is:

[0188] ;

[0189] in, This represents the number of neighboring grid cells; For grid The neighborhood set; For neighborhood grid The variable value.

[0190] In this embodiment of the invention, S2 includes the following sub-steps:

[0191] S21. Based on the preprocessed basic data, determine several growth stages according to the target crop type and identify crop phenology;

[0192] S22. Based on crop phenology, construct physiological regulation functions, specifically including phenological regulation function, radiation regulation function, temperature regulation function, saturated water vapor pressure deficit regulation function, and soil moisture content regulation function;

[0193] S23. Determine stomatal conductance based on phenological regulation function, radiation regulation function, temperature regulation function, saturated water vapor pressure deficit regulation function, and soil moisture content regulation function.

[0194] This step is used to identify the crop growth stage at each grid cell and year, calculate the regulatory parameters related to stomatal opening and physiological activity, and provide a dynamic physiological background for calculating air pollution mechanism quantities.

[0195] Based on the target crop type, several growth stages are pre-defined. These growth stages can be divided according to the growth patterns of different crops. For example, winter wheat can be divided into the pre-wintering stage, the greening and jointing stage, the heading and booting stage, and the grain-filling and ripening stage; corn can be divided into the seedling stage, the jointing stage, the tasseling and silking stage, and the grain-filling and ripening stage. Other crops can be divided into stages according to the same principle. Each grid and each year corresponds to a set of stage start dates and stage end dates.

[0196] Phenological stages are identified by comprehensively considering sowing information, effective accumulated temperature, and remote sensing growth curves. The cumulative effective accumulated temperature is calculated based on the sowing date and crop base temperature; when the cumulative effective accumulated temperature reaches a certain stage threshold, the corresponding growth stage is determined.

[0197] In this embodiment of the invention, in S21, the method for determining the growth stage is as follows: The cumulative effective accumulated temperature is calculated based on the sowing date and the crop base temperature; when the cumulative effective accumulated temperature reaches a threshold, it is determined that the plant has entered the growth stage; the cumulative effective accumulated temperature... The expression is:

[0198] ;

[0199] ;

[0200] in, The average daily temperature; Base temperature for crops; The date sequence number for accumulated temperature; For grid ,years Next stage The starting date sequence number; For grid ,years Next stage End date sequence number; For the stage The corresponding set of dates;

[0201] when Reaching the stage threshold At that time, it is determined that the stage has been entered. . For the stage The accumulated temperature threshold, ℃·d. To improve the accuracy of regional scale identification, the start and end times of the stage can be corrected by combining the time series variation characteristics of NDVI, LAI, or SIF.

[0202] Through the above processing, the results of the growth stage division of the target crop under each grid and each year, as well as the date boundaries of each stage, are obtained, which serve as the basis for subsequent physiological regulation function calculations.

[0203] Based on crop phenology, radiation, temperature, air aridity, and soil moisture conditions, various physiological regulation functions are calculated, and stomatal conductance is calculated accordingly.

[0204] In S22, the phenological regulation function The expression is:

[0205] ;

[0206] in, This refers to the serial number of the start date of the reproductive period; The date number corresponding to the peak physiological activity; This refers to the serial number of the end date of the reproductive period;

[0207] When the crop is in a stage of high physiological activity, the phenological regulation function takes a higher value; when it is in a stage of low physiological activity, the function takes a lower value.

[0208] In S22, the radiation regulation function is used to characterize the promoting effect of light conditions on stomatal opening and leaf physiological activity, as well as the process by which this promoting effect gradually approaches saturation as light intensity increases.

[0209] Radiation regulation function The expression is:

[0210] ;

[0211] in, Photosynthetically active radiation; The light intensity is half-saturated.

[0212] As photosynthetically active radiation increases, the radiation regulation function gradually increases; when radiation reaches a high level, its increase slows down.

[0213] In S22, the temperature regulation function is used to characterize the regulatory effect of temperature conditions on crop physiological activities and stomatal exchange capacity.

[0214] Temperature regulation function The expression is:

[0215] ;

[0216] in, The average daily temperature; The minimum suitable temperature for crop growth; The optimal temperature for crops; This is the optimal temperature for crop growth.

[0217] When the temperature is close to the optimal temperature, the temperature regulation function takes a higher value; when the temperature deviates from the optimal range, the function decreases.

[0218] In S22, the saturated vapor pressure deficit (VPD) regulation function is used to characterize the stomatal closure process that occurs in crops to reduce transpiration water loss when air dryness increases.

[0219] Saturated vapor pressure deficit adjustment function The expression is:

[0220] ;

[0221] in, This indicates a saturated vapor pressure deficit. The threshold for low dryness is [not specified]. The high dryness threshold;

[0222] The function takes a higher value when the air is humid; the function gradually decreases as the air becomes drier.

[0223] In S22, the soil moisture content regulation function is used to characterize the constraint effect of root zone water supply conditions on leaf stomatal opening and crop physiological activity.

[0224] Soil moisture content adjustment function The expression is:

[0225] ;

[0226] in, This refers to the soil moisture content in the root zone. This refers to the moisture content at the wilting point. Field holding capacity;

[0227] In S23, stomatal conductance is used to comprehensively characterize the strength of the leaf gas exchange channels of the target crop under the current phenological stage and environmental conditions. It is a prerequisite parameter for calculating the effective absorbed dose of gaseous pollutants such as ozone, sulfur dioxide, and nitrogen dioxide.

[0228] Pore ​​conductance The expression is:

[0229] ;

[0230] in, This represents the maximum porosity.

[0231] In this embodiment of the invention, S3 includes the following sub-steps:

[0232] S31. Based on stomatal conductance, cumulative effect of ozone concentration exceeding the plant sensitivity threshold, weighted cumulative effect of high-concentration ozone, and effective ozone dose exceeding the flux threshold, the comprehensive ozone mechanism index is determined by weighting.

[0233] S32. Calculate the ozone photosynthesis inhibition index based on the comprehensive ozone mechanism index and the maximum stomatal conductance.

[0234] S33. Calculate the volume extinction coefficient of particulate matter based on the concentrations of PM2.5 and PM10.

[0235] S34. Calculate direct radiation based on the volume extinction coefficient of particulate matter;

[0236] S35. Calculate the scattered radiation;

[0237] S36. Calculate the particulate matter radiation regulation index based on direct radiation and scattered radiation;

[0238] S37. Calculate the particulate matter settling amount;

[0239] S38. Calculate the amount of sulfur dioxide absorbed by the leaves and the amount of sulfur deposition.

[0240] S39. Calculate the soil pH at the end of the stage based on the amount of sulfur deposition, and determine the soil acidification index.

[0241] S311. Calculate the leaf absorption dose of nitrogen dioxide and the amount of nitrogen deposition;

[0242] S312. The comprehensive ozone mechanism index, ozone photosynthesis inhibition index, particulate matter radiation regulation index, and particulate matter deposition amount are used as air quality mechanism indexes, and the leaf absorption dose of sulfur dioxide, sulfur deposition amount, soil acidification index, leaf absorption dose of nitrogen dioxide, and nitrogen deposition amount are used as pollutant mechanism indexes.

[0243] Based on the physiological regulation function and stomatal conductance obtained in step S2, and combined with air pollutant concentration data, mechanistic indicators of the effects of ozone, particulate matter, sulfur dioxide, and nitrogen dioxide on crop growth are constructed to realize the conversion from external pollution levels to the actual intensity of the effects on crops. For gaseous pollutants such as ozone, sulfur dioxide, and nitrogen dioxide, their exposure levels, leaf absorbed doses, or effective fluxes are calculated; for particulate matter, their radiation modulation effect and leaf deposition effect are calculated. After stage aggregation, the mechanistic indicators of various pollutants are used as input features for the crop yield prediction model in step S4.

[0244] Calculate the ozone exposure intensity, stomatal uptake intensity, and combined stress effects on crops. Ozone has a dual effect on crops: on the one hand, the persistent presence of high concentrations of ozone during the growth stage leads to cumulative exposure; on the other hand, ozone must enter the plant through leaf stomata, causing oxidative stress and inhibition of photosynthesis. Therefore, this step considers both exposure and flux indicators.

[0245] In this embodiment of the invention, S31 refers to the cumulative effect of ozone concentration exceeding the plant sensitivity threshold. The expression is:

[0246] ;

[0247] in, This refers to the hourly ozone concentration. For grid ,years ,stage Daytime hourly gatherings;

[0248] In S31, the weighted cumulative effect of high concentrations of ozone The expression is:

[0249] ;

[0250] in, It is an exponential function;

[0251] In S31, the effective ozone dose after exceeding the flux threshold The expression is:

[0252] ;

[0253] in, Unit conversion factor; This refers to the ozone flux threshold. The step size is in hours; Pore ​​conductance;

[0254] In S31, the comprehensive ozone mechanism index The expression is:

[0255] ;

[0256] in, This is the standardized AOT40; This is the standardized W126; For standardization ; The weights of the standardized AOT40; W126 is the weighted and standardized version; For standardization The weights;

[0257] In S32, the ozone photosynthesis inhibition index The expression is:

[0258] ;

[0259] in, Maximum porosity;

[0260] In S33, the volume extinction coefficient of particulate matter is calculated based on the concentrations of PM2.5 and PM10.

[0261] Volume extinction coefficient of particulate matter The expression is:

[0262] ;

[0263] in, PM2.5 extinction efficiency; PM10 extinction efficiency; PM2.5 concentration; PM10 concentration;

[0264] In S34, the direct radiation reaching the canopy under the action of particulate matter is calculated based on the particulate matter extinction coefficient and potential radiation.

[0265] direct radiation The expression is:

[0266] ;

[0267] in, Potential radiation; Equivalent atmospheric optical path;

[0268] In S35, scattered radiation The expression is:

[0269] ;

[0270] in, The particulate scattering coefficient;

[0271] In S36, a particulate matter radiation modulation index is constructed based on the stage-average scattered radiation, stage-average direct radiation, and stage-average potential radiation.

[0272] Particulate matter radiation regulation index The expression is:

[0273] ;

[0274] in, The average scattered radiation of the stage; The average potential radiation for the stage; The average direct radiation during the phase;

[0275] In S37, the particulate matter settling flux is accumulated over time during the growth stage to obtain the stage particulate matter settling amount.

[0276] Particulate matter sedimentation The expression is:

[0277] ;

[0278] ;

[0279] in, The daily time step; This refers to particulate matter settling flux; PM2.5 deposition rate; PM10 deposition velocity; For the stage The corresponding set of dates;

[0280] In S38, the cumulative dose of sulfur dioxide absorbed by the leaves during the stages is calculated based on the concentration of sulfur dioxide in the air, stomatal conductance, and exposure time.

[0281] Sulfur dioxide leaf absorption dose The expression is:

[0282] ;

[0283] in, This refers to the concentration of sulfur dioxide.

[0284] In S38, sulfur deposition The expression is:

[0285] ;

[0286] in, This refers to the settling velocity of sulfur dioxide. This refers to the concentration of sulfur dioxide. For the stage The corresponding set of dates;

[0287] In S39, the soil pH at the end of the stage is updated based on the amount of sulfur deposition in each stage.

[0288] soil pH at the end of the stage The expression is:

[0289] ;

[0290] in, This represents the soil pH at the end of the previous stage. The corresponding coefficient for acidification;

[0291] In S39, soil acidification index The expression is:

[0292] ;

[0293] in, For reference soil pH;

[0294] In S311, nitrogen deposition The expression is:

[0295] ;

[0296] in, This refers to the settling velocity of nitrogen dioxide. This refers to the concentration of nitrogen dioxide.

[0297] In S311, the cumulative dose of nitrogen dioxide absorbed by the leaves during the stages is calculated based on the concentration of nitrogen dioxide in the air, stomatal conductance, and exposure time.

[0298] Nitrogen dioxide leaf absorption dose The expression is:

[0299] ;

[0300] in, This represents the concentration of nitrogen dioxide.

[0301] In this embodiment of the invention, S4 includes the following sub-steps:

[0302] S41. Construct feature vectors based on air quality mechanism indicators and pollutant mechanism indicators;

[0303] S42. Based on the feature vectors, construct a crop yield prediction model and determine the predicted crop yield;

[0304] S43. Supervised learning training of crop yield prediction models using historical samples;

[0305] S44. Set the air quality mechanism index to the pollution-free baseline level to obtain the baseline output;

[0306] S45. Calculate the total loss based on the crop yield prediction model trained under supervised learning and the baseline yield;

[0307] S46. Calculate the relative loss rate based on the total loss.

[0308] S47. Calculate the pollutant contribution rate;

[0309] S46. Calculate the total spatial aggregate yield based on the crop yield forecast;

[0310] S47. Output the predicted crop yield, baseline yield, total loss, relative loss rate, pollutant contribution rate, and spatial aggregate total yield as results.

[0311] The air quality mechanism indicators obtained in step S3 are combined with meteorological, remote sensing, soil and management features to construct crop yield prediction input features. Machine learning models are used to obtain raster-scale crop yield prediction results, and the yield loss, relative loss rate, pollutant contribution rate and administrative unit scale aggregated results caused by air pollution are calculated.

[0312] In this embodiment of the invention, in S41, the feature vector The expression is:

[0313] ;

[0314] in, This represents the cumulative effect of ozone concentration exceeding the plant sensitivity threshold. This represents the weighted cumulative effect of high concentrations of ozone. The effective ozone dose after exceeding the flux threshold; As a comprehensive indicator of ozone mechanism; As an indicator of ozone inhibition of photosynthesis; As an indicator of particulate matter radiation regulation; This refers to particulate matter settling volume; This refers to the leaf absorption dose of sulfur dioxide; This refers to the amount of sulfur deposition. Indicators of soil acidification; This refers to the leaf absorption dose of sulfur dioxide; This refers to nitrogen deposition. The average temperature of the stage; This refers to the cumulative precipitation over a given period. The average photosynthetically active radiation (APRA) during the phase; This represents the average soil moisture content for that stage. The leaf area index is the stage index; For stage NDVI; For staged solar-induced chlorophyll fluorescence; This is a soil attribute vector; To manage feature vectors;

[0315] In step S42, a machine learning crop yield prediction model is established based on the feature vector constructed in step S41 to obtain the raster-scale predicted yield per unit area.

[0316] Crop yield forecast The expression is:

[0317] ;

[0318] in, For machine learning model functions;

[0319] The machine learning model function can be one or more of the following: random forest, XGBoost, LightGBM, support vector regression, deep neural network, or ensemble model.

[0320] In S43, the expression for the objective function for supervised learning training is:

[0321] ;

[0322] in, For training the objective function; For the set of model parameters; This represents the actual crop yield; The regularization coefficient is used. It is a regularization term;

[0323] During model training, cross-validation, spatiotemporal hierarchical validation, or hyperparameter optimization can be combined to improve the model's generalization ability. Specifically: .

[0324] In S44, the baseline production The expression is:

[0325] ;

[0326] in, These are the feature vectors under no-pollution or baseline-pollution scenarios;

[0327] In S45, the total loss The expression is:

[0328] ;

[0329] In S46, the relative loss level caused by air pollution is calculated to characterize the relative impact of air pollution on crop yield, and can be used for subsequent risk classification and spatial comparison analysis.

[0330] relative loss rate The expression is:

[0331] ;

[0332] In step S47, the contribution of different pollutants to crop yield changes or losses is calculated. Machine learning interpretation methods, such as SHAP, permutation importance, or feature attribution methods, are preferably used to obtain the contribution values ​​of each pollutant's related features to the yield prediction results.

[0333] Pollutant contribution rate The expression is:

[0334] ;

[0335] in, Index for pollutant categories; pollutants The absolute contribution value in production forecasting;

[0336] In S46, the raster-scale yield results are aggregated to administrative units or other statistical units to obtain the total yield results.

[0337] Spatial Aggregation Total Output The expression is:

[0338] ;

[0339] in, For target crops in the grid ,years The planting area.

[0340] Calculate the area-weighted average loss rate at the administrative unit level. The average loss rate per administrative unit is calculated using the following formula:

[0341] ;

[0342] in, For administrative units The average loss rate, %.

[0343] Those skilled in the art will recognize that the embodiments described herein are intended to help the reader understand the principles of the invention, and should be understood that the scope of protection of the invention is not limited to such specific statements and embodiments. Those skilled in the art can make various other specific modifications and combinations based on the technical teachings disclosed in this invention without departing from the spirit of the invention, and these modifications and combinations are still within the scope of protection of this invention.

Claims

1. A grid-scale method for predicting crop yield based on the mechanism of air quality effects, characterized in that, Includes the following steps: S1. Divide the study area into grids and collect basic data. Then, preprocess the basic data in sequence. S2. Based on the preprocessed basic data, construct the physiological regulation function and determine the stomatal conductance; S3. Determine air quality mechanism indicators and pollutant mechanism indicators based on porosity; S4. Based on air quality mechanism indicators and pollutant mechanism indicators, determine the predicted crop yield and output the baseline yield, total loss, relative loss rate, pollutant contribution rate, and spatial aggregate total yield.

2. The method for predicting crop yield based on the air quality mechanism at the grid scale according to claim 1, characterized in that, In S1, the preprocessing includes time scale unification, spatial scale unification, and missing value imputation; The expression for the unified time scale is: ； in, This represents the daily average value. These are the original variable values ​​on an hourly scale; The number of valid hours for the day; The expression for the unified spatial scale is: ； ； in, These are the interpolated raster variable values; For the site Observed values; For grid With the site The weights; For grid With the site The distance; This is the distance attenuation coefficient; The total number of sites participating in the interpolation; The missing value filling specifically refers to: If there are missing dates or hours, the expression for filling them in is: ； in, The corrected variable values; The variable value is from the previous day; The variable value for the following day; If consecutive missing values ​​exceed a preset threshold, the expression for filling them in is: ； in, This represents the number of neighboring grid cells; For grid The neighborhood set; For neighborhood grid The variable value.

3. The method for predicting crop yield based on the air quality mechanism at the grid scale according to claim 1, characterized in that, S2 includes the following sub-steps: S21. Based on the preprocessed basic data, determine several growth stages according to the target crop type and identify crop phenology; S22. Based on crop phenology, construct physiological regulation functions, specifically including phenological regulation function, radiation regulation function, temperature regulation function, saturated water vapor pressure deficit regulation function, and soil moisture content regulation function; S23. Determine stomatal conductance based on phenological regulation function, radiation regulation function, temperature regulation function, saturated water vapor pressure deficit regulation function, and soil moisture content regulation function.

4. The method for predicting crop yield based on the air quality mechanism at the grid scale according to claim 3, characterized in that, In step S21, the method for determining the growth stage is as follows: The cumulative effective accumulated temperature is calculated based on the sowing date and the crop base temperature; when the cumulative effective accumulated temperature reaches a threshold, it is determined that the crop has entered the growth stage; the cumulative effective accumulated temperature... The expression is: ； ； in, The average daily temperature; Base temperature for crops; The date sequence number for accumulated temperature; For grid ,years Next stage The starting date sequence number; For grid ,years Next stage End date sequence number; For the stage The corresponding set of dates; In S22, the phenological adjustment function The expression is: ； in, This refers to the serial number of the start date of the reproductive period; The date number corresponding to the peak physiological activity; This refers to the serial number of the end date of the reproductive period; In S22, the radiation adjustment function The expression is: ； in, Photosynthetically active radiation; The light intensity is half-saturated. In S22, the temperature adjustment function The expression is: ； in, The average daily temperature; The minimum suitable temperature for crop growth; The optimal temperature for crops; This is the optimal temperature for crop growth. In S22, the saturated water vapor pressure deficit adjustment function The expression is: ； in, This indicates a saturated vapor pressure deficit. Low dryness threshold; The threshold for high dryness; In S22, the soil moisture content adjustment function The expression is: ； in, This refers to the soil moisture content in the root zone. This refers to the moisture content at the wilting point. Field holding capacity; In S23, the porosity The expression is: ； in, This represents the maximum porosity.

5. The method for predicting crop yield based on the air quality mechanism at the grid scale according to claim 1, characterized in that, S3 includes the following sub-steps: S31. Based on stomatal conductance, cumulative effect of ozone concentration exceeding the plant sensitivity threshold, weighted cumulative effect of high-concentration ozone, and effective ozone dose exceeding the flux threshold, the comprehensive ozone mechanism index is determined by weighting. S32. Calculate the ozone photosynthesis inhibition index based on the comprehensive ozone mechanism index and the maximum stomatal conductance. S33. Calculate the volume extinction coefficient of particulate matter based on the concentrations of PM2.5 and PM10. S34. Calculate direct radiation based on the volume extinction coefficient of particulate matter; S35. Calculate the scattered radiation; S36. Calculate the particulate matter radiation regulation index based on direct radiation and scattered radiation; S37. Calculate the particulate matter settling amount; S38. Calculate the amount of sulfur dioxide absorbed by the leaves and the amount of sulfur deposition. S39. Calculate the soil pH at the end of the stage based on the amount of sulfur deposition, and determine the soil acidification index. S311. Calculate the leaf absorption dose of nitrogen dioxide and the amount of nitrogen deposition; S312. The comprehensive ozone mechanism index, ozone photosynthesis inhibition index, particulate matter radiation regulation index, and particulate matter deposition amount are used as air quality mechanism indexes, and the leaf absorption dose of sulfur dioxide, sulfur deposition amount, soil acidification index, leaf absorption dose of nitrogen dioxide, and nitrogen deposition amount are used as pollutant mechanism indexes.

6. The grid-scale crop yield prediction method based on air quality mechanisms according to claim 5, characterized in that, In S31, the cumulative effect of ozone concentration exceeding the plant sensitivity threshold The expression is: ； in, This refers to the hourly ozone concentration. For grid ,years ,stage Daytime hourly gatherings; In S31, the weighted cumulative effect of high concentration ozone The expression is: ； in, It is an exponential function; In S31, the effective ozone dose after exceeding the flux threshold The expression is: ； in, Unit conversion factor; This refers to the ozone flux threshold. The step size is in hours; Pore ​​conductance; In S31, the comprehensive ozone mechanism index The expression is: ； in, This is the standardized AOT40; This is the standardized W126; For standardization ; The weights of the standardized AOT40; W126 is the weighted standardization; For standardization The weights; In S32, the ozone photosynthesis inhibition index The expression is: ； in, Maximum porosity; In S33, the particulate matter volume extinction coefficient The expression is: ； in, PM2.5 extinction efficiency; PM10 extinction efficiency; PM2.5 concentration; PM10 concentration; In S34, direct radiation The expression is: ； in, Potential radiation; Equivalent atmospheric optical path; In S35, scattered radiation The expression is: ； in, The particulate scattering coefficient; In S36, the particulate matter radiation regulation index The expression is: ； in, The average scattered radiation of the stage; The average potential radiation for the stage; The average direct radiation during the phase; In S37, the particulate matter sedimentation amount The expression is: ； ； in, The daily time step; This refers to particulate matter settling flux; PM2.5 deposition rate; PM10 deposition velocity; For the stage The corresponding set of dates; In S38, the sulfur dioxide leaf absorption dose The expression is: ； in, This refers to the concentration of sulfur dioxide. In S38, the amount of sulfur deposition The expression is: ； in, This refers to the settling velocity of sulfur dioxide. This refers to the concentration of sulfur dioxide. For the stage The corresponding set of dates; In S39, the soil pH at the end of the stage The expression is: ； in, This represents the soil pH at the end of the previous stage. The corresponding coefficient for acidification; In S39, the soil acidification index The expression is: ； in, For reference soil pH; In S311, the nitrogen deposition amount The expression is: ； in, This refers to the settling velocity of nitrogen dioxide. This refers to the concentration of nitrogen dioxide. In S311, the nitrogen dioxide leaf absorption dose The expression is: ； in, This represents the concentration of nitrogen dioxide.

7. The method for predicting crop yield based on the air quality mechanism at the grid scale according to claim 1, characterized in that, S4 includes the following sub-steps: S41. Construct feature vectors based on air quality mechanism indicators and pollutant mechanism indicators; S42. Based on the feature vectors, construct a crop yield prediction model and determine the predicted crop yield; S43. Supervised learning training of crop yield prediction models using historical samples; S44. Set the air quality mechanism index to the pollution-free baseline level to obtain the baseline output; S45. Calculate the total loss based on the crop yield prediction model trained under supervised learning and the baseline yield; S46. Calculate the relative loss rate based on the total loss. S47. Calculate the pollutant contribution rate; S46. Calculate the total spatial aggregate yield based on the crop yield forecast; S47. Output the predicted crop yield, baseline yield, total loss, relative loss rate, pollutant contribution rate, and spatial aggregate total yield.

8. The method for predicting crop yield based on the air quality mechanism at the grid scale according to claim 7, characterized in that, In S41, the feature vector The expression is: ； in, This represents the cumulative effect of ozone concentration exceeding the plant sensitivity threshold. This represents the weighted cumulative effect of high concentrations of ozone. The effective ozone dose after exceeding the flux threshold; As a comprehensive indicator of ozone mechanism; As an indicator of ozone inhibition of photosynthesis; As an indicator of particulate matter radiation regulation; This refers to particulate matter settling volume; This refers to the leaf absorption dose of sulfur dioxide; This refers to the amount of sulfur deposition. Indicators of soil acidification; This refers to the leaf absorption dose of sulfur dioxide; This refers to nitrogen deposition. The average temperature of the stage; This refers to the cumulative precipitation over a given period. The average photosynthetically active radiation (APRA) during the phase; This represents the average soil moisture content for that stage. The leaf area index is the stage index; For stage NDVI; For staged solar-induced chlorophyll fluorescence; This is a soil attribute vector; To manage feature vectors; In S42, crop yield prediction The expression is: ； in, For machine learning model functions; In step S43, the expression for the objective function for supervised learning training is: ； in, For training the objective function; For the set of model parameters; This represents the actual crop yield; The regularization coefficient is used. It is a regularization term; In S44, the benchmark output The expression is: ； in, These are the feature vectors under no-pollution or baseline-pollution scenarios; In S45, the total loss amount The expression is: ； In S46, the relative loss rate The expression is: ； In S47, the pollutant contribution rate The expression is: ； in, Index for pollutant categories; pollutants The absolute contribution value in production forecasting; In S46, the total spatial aggregation yield The expression is: ； in, For target crops in the grid ,years The planting area.

Images