Method for evaluating effects of surface ozone on winter wheat yield based on multi-factor stress

Abstract

The present application relates to the technical field of crop state evaluation, and particularly relates to a method for evaluating the influence of ground ozone on winter wheat yield based on multi-factor stress, comprising: dividing a research area to be evaluated according to the positions of crop meteorological stations in the research area to be evaluated, and obtaining each sub-research area in the research area to be evaluated; obtaining a basic data file of each sub-research area, and determining an optimal model parameter group corresponding to an O3modified-NWheat model; and determining potential yield loss caused by ozone to each sub-research area according to the O3modified-NWheat model corresponding to the optimal model parameter group and the basic data file of each sub-research area. The present application can realize accurate evaluation of potential yield loss of each sub-research area caused by ozone, and can help to provide a theoretical basis for formulating a regional management scheme for farmland ozone in a winter wheat main production area.

Description

Technical Field

[0001] This invention relates to the field of crop condition assessment technology, specifically to a method for assessing the impact of surface ozone on winter wheat yield based on multi-factor stress. Background Technology

[0002] Ozone is generated through photochemical reactions from volatile organic compounds, nitrogen oxides, carbon monoxide, and methane emitted by biological sources and human activities. It is a secondary atmospheric pollutant with significant negative impacts on surface vegetation. Investigating the impact of ozone on the yield of major food crops is crucial for ensuring food security. Wheat is significantly more sensitive to ozone than other food crops. Near-surface ozone can ultimately reduce wheat yield by affecting its photosynthetic rate and accelerating leaf senescence. Winter wheat producing areas, as major grain production bases, account for more than three-quarters of the total annual wheat planting area and yield. Therefore, these areas are the most severely affected by ozone pollution. Assessing the extent of ozone damage to winter wheat production yield is necessary to provide a theoretical basis for setting subsequent targets for ozone pollution control in farmland.

[0003] Existing methods quantify the relationship between wheat yield and ozone concentration by setting different ozone dosage indicators. Among the numerous ozone assessment indicators, M7 (daily average ozone concentration from 09:00 to 16:00), AOT40 (cumulative value of hourly ozone concentration exceeding 40 ppb during the day), and POD are used. Y Indicators such as ozone uptake flux are widely used in the current regional assessment of the impact of ozone on wheat yield. However, the use of ozone assessment indicators is often limited when assessing the impact of ozone on wheat yield at a regional scale. More comprehensive methods and approaches for assessing the damage of regional ozone pollution to crops still need to be explored. For example, ozone indicators often cannot fully consider the combined effects of ozone and other abiotic stresses; wheat planting environments and field management practices vary in different regions, and using a single ozone assessment indicator to assess the degree of ozone impact on winter wheat has poor robustness, increasing the uncertainty of regional winter wheat status assessment.

[0004] Furthermore, existing crop models that comprehensively simulate the relationship between crop growth and development and crop variety, environmental conditions, and field management practices are widely used in studies simulating and analyzing the impact of changes in crop growth environmental factors on crop growth and yield. Among these, the Decision Support System for Agro-technology Transfer (DSSAT) and the Agricultural Production Systems Simulator (APSIM) are currently the most widely used wheat crop models. After integrating the wheat simulation module from APSIM into the DSSAT model, the DSSAT model, which includes three wheat sub-models (CERES, CROPSIM, and NWheat), demonstrates superior performance in simulating wheat compared to other crop models. However, although many studies have confirmed the negative impact of ozone on wheat yield, the lack of incorporation of the mechanism of ozone's influence on wheat state into the model, incomplete basic data files, and failure to localize model training parameters according to specific study areas result in low accuracy of model outputs, meaning that the extent of ozone's potential yield loss cannot be accurately assessed. Summary of the Invention

[0005] To address the aforementioned technical problem of the inability to accurately assess the extent of ozone damage to potential crop yields, the present invention aims to provide a method for assessing the impact of surface ozone on winter wheat yield based on multi-factor stress. The specific technical solution adopted is as follows: One embodiment of the present invention provides a method for assessing the impact of surface ozone on winter wheat yield based on multi-factor stress, the method comprising the following steps: Based on the location of crop meteorological stations within the study area to be evaluated, the study area to be evaluated is divided into sub-study areas within the study area to be evaluated. Obtain the basic data files for each sub-study area and determine the optimal model parameter set corresponding to the O3 modified-NWheat model; Based on the O3 modified-NWheat model corresponding to the optimal model parameter set and the basic data files of each sub-study area, the potential yield loss caused by ozone in each sub-study area is determined.

[0006] Furthermore, the basic data files include variety files, soil files, meteorological files, and field management files. The variety files are winter wheat varieties from various sub-study areas during the study period. The soil files include station latitude and longitude, soil type, gravel content, soil color, surface slope, organic carbon content, pH value, soil moisture content, saturated moisture content, soil saturated hydraulic conductivity, field water holding capacity, and wilting point water holding capacity. The meteorological files include daily maximum temperature, daily minimum temperature, precipitation, solar radiation, wind speed, relative humidity, CO2 concentration, evapotranspiration, and target ozone index. The field management files include management data and growth and yield record data. The management data includes land preparation method, sowing date, sowing method, sowing depth, sowing density, irrigation date, irrigation amount, fertilization date, fertilizer type, and fertilizer amount. The growth and yield record data includes winter wheat phenological indicators, biomass, grain weight, thousand-grain weight, and actual yield.

[0007] Furthermore, the optimal model parameter set corresponding to the O3 modified-NWheat model is determined, including: Based on the input data required for the O3 modified-NWheat model, create a first preset number of variety genetic parameters in the variety genetic parameter file; wherein, a set of variety genetic parameters includes a second preset number of crop genetic parameters and a third preset number of soil phenotypic parameters; The GLUE method was used to optimize the genetic parameters of each group of varieties until the GLUE was run a fourth preset number of times, and the model output data corresponding to the optimized genetic parameters of each group of varieties was obtained; wherein, the model output data is the simulated value of crop status in each sub-study area; Based on the simulated and observed crop state values ​​corresponding to the genetic parameters of each variety, the weights, mean, and variances of the genetic parameters of each variety are calculated. Then, based on the weights, mean, and variances, the optimal model parameter set corresponding to the O3 modified-NWheat model is determined.

[0008] Furthermore, the process of calculating the weights, mean, and variance of the genetic parameters for each group of varieties includes: Based on the simulated and observed crop state values ​​corresponding to the genetic parameters of each variety, calculate the likelihood value between the simulated and observed crop state values. Based on the likelihood values ​​corresponding to the genetic parameters of each group of varieties, calculate the weight, mean, and variance of the genetic parameters of each group of varieties.

[0009] Furthermore, the formula for calculating the weights of the genetic parameters for each group of varieties is as follows: in, Let be the weights of the genetic parameters of the i-th variety group. Let N be the likelihood value corresponding to the genetic parameters of the i-th variety, N be the number of the genetic parameter groups of the variety, and i be the sequence number of the genetic parameter group of the variety. The formula for calculating the mean of genetic parameters for each group of varieties is: in, The mean values ​​of the genetic parameters for each group of varieties are given. The weights of the likelihood values ​​corresponding to the genetic parameters of the i-th variety group are: Let N be the genetic parameters of the i-th variety group, N be the number of the variety genetic parameter groups, and i be the sequence number of the variety genetic parameter group. The formula for calculating the variance of genetic parameters for each group of varieties is: in, The variance of the genetic parameters for each group of varieties is given. The weights of the likelihood values ​​corresponding to the genetic parameters of the i-th variety group are: Let N be the genetic parameters of the i-th variety group, and N be the number of variety genetic parameter groups. This represents the mean of the genetic parameters for each group of varieties.

[0010] Furthermore, based on the location of crop meteorological stations within the study area to be evaluated, the study area is divided into sub-study areas, including: Using the Thiessen polygon method, the study area to be evaluated is divided into sub-study areas based on the location of crop meteorological stations within the study area to be evaluated. Among them, the edge distance of the sub-study area is the shortest straight-line distance to the corresponding crop meteorological station.

[0011] Furthermore, based on the optimal model parameter set corresponding to the O3 modified-NWheat model and the basic data files of each sub-study area, the potential yield loss caused by ozone in each sub-study area was determined, including: The basic data files of each sub-study area are used as input data and input into the O3 modified-NWheat model corresponding to the optimal model parameter set to determine the crop yield of each sub-study area under the influence of ozone. The O3 modified-NWheat model is then used to simulate the crop yield when it is not affected by ozone. Based on the crop yield of each sub-study area under the influence of ozone and the crop yield when it is not affected by ozone, the difference between the crop yield when it is not affected by ozone and the crop yield when it is affected by ozone is determined as the potential yield loss caused by the corresponding sub-study area.

[0012] The present invention has the following beneficial effects: This invention provides a method for assessing the impact of surface ozone on winter wheat yield based on multi-factor stress. First, the NWheat model is chosen as an improvement premise for the O3 modified-NWheat model because ozone has a significant impact on wheat leaf senescence. The uncertainty brought about by the model's architecture mainly stems from its consideration of crop growth mechanisms. For example, the CERES and CROPSIM models do not consider the senescence process of wheat leaves. The uncertainty caused by the model structure is often greater than the uncertainty caused by the model parameters. Therefore, choosing the NWheat model is beneficial to improving the accuracy of the model output results. Second, the localization of model parameters is crucial for the successful application of the model in different locations. Existing research shows that model parameter localization can be achieved by selecting more sensitive parameters based on sensitivity analysis. However, selecting the appropriate number of parameters and which sensitive parameters to calibrate remains a challenge in model parameter localization. Based on existing research, this invention selects eight varietal genetic parameters and three soil parameters that are relatively sensitive to wheat growth period and yield for localization, enabling the model to achieve excellent simulation performance. This verifies the rationality of the selected parameters to a certain extent and improves the accuracy of the model output results. Furthermore, the availability and errors of model input data also contribute to the uncertainty of model simulation results. High-quality spatial input data that can well represent crop growth conditions in the study area can significantly reduce this uncertainty. Therefore, on the one hand, more environmental factors were added to the minimum dataset required for model operation, such as wind speed and near-surface relative humidity, to improve the model's accuracy in simulating surface evapotranspiration. On the other hand, the input file processing was optimized to minimize random errors caused by model operators. Finally, considering the three main sources of model uncertainty, the GLUE parameter optimization method was used to localize the model parameters. The ozone-modified NWheat model was successfully used to simulate the potential yield loss caused by ozone in the main winter wheat producing areas of my country during the study period. This invention accurately assesses the degree of potential yield loss caused by ozone, provides a theoretical basis for formulating regional farmland ozone pollution prevention and control measures, and can serve to improve the production potential of winter wheat in the study area. Attached Figure Description

[0013] To more clearly illustrate the technical solutions and advantages in the embodiments of the present invention or the prior art, the drawings used in the description of the embodiments or the prior art will be briefly introduced below. Obviously, the drawings described below are only some embodiments of the present invention. For those skilled in the art, other drawings can be obtained based on these drawings without creative effort.

[0014] Figure 1 This is a flowchart of a method for assessing the impact of surface ozone on winter wheat yield based on multi-factor stress according to the present invention. Figure 2 This is a schematic diagram of the area to be removed from the study area to be evaluated in an embodiment of the present invention; Figure 3 This is a schematic diagram of each sub-study area within the study area to be evaluated in an embodiment of the present invention; Figure 4 A system framework diagram of the basic data files required for the crop state estimation model in this embodiment of the invention; Figure 5 This document describes the changes in winter wheat yield in major winter wheat producing areas from 2015 to 2018, after removing the effects of ozone. Figure 6 This is a statistical result of the potential loss of winter wheat yield due to ozone in 25 regions within the main winter wheat producing areas of this invention. Figure 7 This is the spatial distribution of winter wheat yield recorded by agricultural meteorological stations in 25 sub-regions of the main winter wheat producing areas in this embodiment of the invention. Figure 8 This is a training result diagram of the flowering period in an embodiment of the present invention; Figure 9 This is a diagram showing the verification results of the flowering period in an embodiment of the present invention; Figure 10 This is the spatial distribution of the error between observed and simulated values ​​of winter wheat flowering period in the study area in this embodiment of the invention; Figure 11 This is a training result diagram of the mature stage in an embodiment of the present invention; Figure 12 This is a graph showing the verification results of the maturity stage in an embodiment of the present invention; Figure 13 This is the spatial distribution of the error between observed and simulated values ​​of winter wheat maturity in the study area in this embodiment of the invention; Figure 14 This is a graph showing the training results of the model in this embodiment of the invention for simulating winter wheat yield; Figure 15 This is a graph showing the verification results of the model in this embodiment of the invention for simulating winter wheat yield; Figure 16 This represents the spatial distribution of the error between the observed and simulated values ​​of winter wheat yield in the study area in this embodiment of the invention. Detailed Implementation

[0015] To further illustrate the technical means and effects adopted by the present invention to achieve its intended purpose, the specific implementation methods, structures, features, and effects of the technical solution proposed according to the present invention are described in detail below with reference to the accompanying drawings and preferred embodiments. In the following description, different "one embodiment" or "another embodiment" do not necessarily refer to the same embodiment. Furthermore, specific features, structures, or characteristics in one or more embodiments can be combined in any suitable form.

[0016] Unless otherwise defined, all technical and scientific terms used herein have the same meaning as commonly understood by one of ordinary skill in the art to which this invention pertains.

[0017] Current methods for assessing ozone damage to crop yield largely rely on concentration response functions obtained from representative experiments. However, these methods often fail to capture spatial differences among various environmental factors and, to some extent, ignore the combined effects of multiple abiotic stresses, leading to low accuracy in estimating crop yield. To overcome the shortcomings of these methods in assessing ozone damage to crop yield, which suffer from poor robustness and low accuracy in estimating potential crop yield loss, this embodiment provides a method for assessing the impact of surface ozone on winter wheat yield based on multi-factor stress, such as... Figure 1 As shown, it includes the following steps: S1. Based on the location of crop meteorological stations within the study area to be evaluated, the study area to be evaluated is divided into sub-study areas within the study area to be evaluated.

[0018] In this embodiment, to more rationally zonate varieties and management practices based on existing winter wheat planting areas, a database of wheat varieties and field management information within the study area to be evaluated was established, using 25 representative agricultural meteorological stations. This database provides a reasonable data foundation for the evaluation of the study area. The number of agricultural meteorological stations is highly correlated with the planting areas of different grain crops within the study area, indicating that areas without winter wheat agricultural meteorological stations are not necessarily major winter wheat producing areas. For example, northern regions such as Shanxi and Hebei are mostly mountainous, and while there are a few winter wheat planting areas in the north, the varieties planted and the field management methods differ significantly from those in other regions. Therefore, to improve the rationality of the zoning and reduce the uncertainty brought about by different varieties and field management methods to yield estimation, areas in the north without agricultural meteorological stations and with predominantly mountainous terrain were excluded. A schematic diagram of the excluded areas within the study area to be evaluated is shown below. Figure 2 As shown, the excluded areas include Taiyuan, Yangquan, Xinzhou, Datong, Zhangjiakou, and Chengde.

[0019] After obtaining the excluded study area, the Thiessen polygon method was used. Based on the location of agricultural meteorological stations within the study area to be evaluated, ArcGIS (Geographic Information System software) was used to divide the study area into sub-study areas. The implementation process of the Thiessen polygon method and the GIS software is existing technology and not within the scope of this invention; therefore, it will not be described in detail here. Each sub-study area exhibits the characteristic that any location within the polygon has the shortest straight-line distance to the corresponding agricultural meteorological station. This facilitates the analysis of the distribution of winter wheat variety planting and management measures within the corresponding sub-study area using data collected from agricultural meteorological stations within the sub-study area. A schematic diagram of each sub-study area within the study area to be evaluated is shown below. Figure 3 As shown.

[0020] S2, Obtain the basic data files for each sub-study area and determine the optimal model parameter set corresponding to the O3 modified-NWheat model. The steps include: The first step is to obtain the basic data files for each sub-study area.

[0021] The crop state estimation model (O3 modified-NWheat) mainly comprises three parts: crop model, data file management, and system application analysis. Each part has its own independent function but can also be used in conjunction with others. This model fully considers the combined effects of crop variety genetic characteristics, soil, climate, and management practices on crop growth, and can be used to simulate the impact of different crop planting methods and management practices on crop growth and production. Because the crop state estimation model is based on a semi-empirical, semi-parametric principle, it requires comprehensive crop growth data, most of which must be obtained experimentally or through real-time observation. The basic data files required for the crop state estimation model mainly include crop variety files, soil files, meteorological files, and field management files. The crops can be winter wheat, or other crops such as corn and rice. The system framework diagram of the basic data files required for the crop state estimation model is shown below. Figure 4 As shown. In Figure 4 The documents include crop variety files (.CUL), soil files (.SOL), meteorological files (.WTH), and field management files (.WHX).

[0022] The first sub-step is the crop variety file (.CUL).

[0023] In this embodiment, agricultural meteorological monitoring stations recorded data on winter wheat varieties planted at each station in different years. Some stations used slightly different winter wheat varieties between years. Considering the differences in winter wheat varieties planted in different sub-regions and the feasibility of localizing variety parameters, based on the winter wheat varieties planted from 2015 to 2018, each station selected dominant winter wheat varieties for each sub-study area according to the planting frequency of wheat varieties. If the winter wheat seedling varieties planted from 2015 to 2018 were different, data on varieties planted in previous years were added for reference. It is worth noting that, due to the lack of data on ozone tolerance indicators for each winter wheat variety, to prevent exaggerating the impact of ozone on the state of winter wheat, it was assumed that all winter wheat varieties planted in all sub-study areas were ozone-tolerant varieties.

[0024] The second sub-step is the soil file (.SOL).

[0025] In this embodiment, due to the large amount of soil file information required by the crop state estimation model, mainly including a series of parameters such as station latitude and longitude, soil type, gravel content, soil color, surface slope, organic carbon content, pH value, soil moisture content, and saturated water content, there is a lack of comprehensive soil input data with different resolutions. In recent years, establishing a spatially continuous high-precision soil database globally has become one of the key focuses of related scholars. The most representative example is the global 1 km soil grid database released by the International Soil Reference and Information Center in 2014. Subsequently, the International Institute for Climate and Society and the International Food Policy Research Institute, among other institutions, jointly derived hydraulic properties crucial for simulating crop growth, such as soil saturated hydraulic conductivity, field capacity, wilting point, and saturation, using soil transport functions based on this dataset. They also developed a 10 km resolution soil dataset (Global High-Resolution Soil Profile Database for Crop Modeling Applications) compatible with the DSSAT crop model by combining the HarvestChoice (HC27) soil database. It has been widely used in the field of large-scale crop model yield estimation research, and this dataset also provides the most reliable soil data information for the soil files required in this embodiment.

[0026] The third sub-step is the weather file (.WTH).

[0027] In this embodiment, to comprehensively consider the impact of meteorological factors on winter wheat growth, in addition to daily maximum temperature, daily minimum temperature, precipitation, and solar radiation, wind speed and relative humidity are also added to calculate evapotranspiration, which helps improve the accuracy of evapotranspiration calculation. Simultaneously, CO2 concentration data, derived from the CO2 basic standard data file included in the DSSAT model, is also taken into account. Evapotranspiration is calculated using the FAO-Penman-Monteith formula, which can be automatically calculated by the model. The formula is as follows: in, Evapotranspiration, in mm / day -1 R n Net radiation on the surface of crops, measured in MJ / m². -2 day -1 G represents soil heat flux density, in MJ / m³. -2 day -1 γ is the humidity calculation constant, with units of kPa℃. -1 T represents the average daily surface temperature in °C, and u2 represents the surface wind speed in m / s². -1 e s This is the saturated vapor pressure, in kPa, e a This is the actual vapor pressure, in kPa. The pressure difference is saturated gas pressure, in kPa; Δ is the slope of the vapor pressure curve, in kPa. .

[0028] The steps for obtaining the target ozone concentration index from the meteorological files (.WTH) of each sub-study area include: In this embodiment, surface ozone monitoring data for each sub-study area can be obtained from an environmental monitoring center. Based on the hourly ozone data, the daily M7 ozone concentration for each sub-study area is calculated, and the M7 ozone concentration is determined as the target ozone concentration index. The calculation process for the daily M7 ozone concentration is prior art and is not within the scope of protection of this invention; therefore, it will not be described in detail here.

[0029] It should be noted that the M7 ozone concentration in meteorological documents plays a crucial role in improving crop state estimation models. The M7 ozone concentration is used to measure the impact of ozone on winter wheat growth and yield for two main reasons: first, numerous experiments have demonstrated that the M7 ozone concentration can be used to assess the impact of ozone stress on wheat yield and can be used to model related studies; second, although some studies use threshold-based cumulative ozone exposure indicators to assess the impact of ozone on wheat yield, these indicators ignore the damage caused by low concentrations of ozone, resulting in lower accuracy in wheat yield damage assessments. The M7 ozone concentration not only takes into account high ozone levels but also shows a significant impact on winter wheat state in the early stages of growth and development. Based on prior knowledge, the ozone stress threshold is set at 25 ppb. When the ozone concentration exceeds the set threshold, photosynthesis will decrease by 0-1 times, and leaves will also accelerate senescence at a rate greater than 1.

[0030] The steps for determining the ozone stress function based on the target ozone concentration index and ozone-related data include: It should be noted that the ozone stress function includes the photosynthetic reduction factor caused by ozone stress without interaction with water and atmospheric CO2 (carbon dioxide) stress, the photosynthetic reduction factor when ozone stress interacts with water and atmospheric CO2 stress, the amount of dry matter accumulation reduced due to decreased photosynthesis caused by growth stress, and the leaf senescence acceleration factor caused by ozone stress; ozone-related data include wheat's tolerance to ozone, the effect of atmospheric CO2 concentration greater than 350 ppm on wheat's potential dry matter, the effect of water stress on photosynthesis, nitrogen stress factor, and temperature stress factor.

[0031] The formula for calculating the photosynthetic decrease factor caused by ozone stress without interaction with moisture and atmospheric CO2 stress is as follows: Among them, FO3 is the photosynthetic reduction factor caused by ozone stress without interaction with moisture and CO2 stress in the air; FOZ1 and FOZ2 are values ​​of 0.0006 and 1.015 (tolerant type), 0.001 and 1.025 (moderately sensitive type), and 0.005 and 1.125 (sensitive type) respectively, based on the ozone tolerance of wheat; and OZONX is the target ozone concentration index.

[0032] The formula for calculating the photosynthetic decrease factor when ozone stress interacts with moisture and atmospheric CO2 stress is as follows: Among them, PRFO3 is the photosynthetic reduction factor under the interaction of ozone stress with moisture and atmospheric CO2 stress, Ruefactor is the effect of atmospheric CO2 concentration greater than 350 ppm on the potential dry matter of wheat, and Swdef(photo nw ) represents the effect of water stress on photosynthesis, while FO3 represents the decrease in photosynthesis caused by ozone stress without interaction with water and atmospheric CO2 stress. To find the minimum value of the function.

[0033] The formula for calculating the decrease in dry matter accumulation due to reduced photosynthesis caused by growth stress is as follows: Where Optfr represents the amount of dry matter accumulation reduced due to decreased photosynthesis caused by growth stress, and Swdef(photo nw The graph represents the effect of water stress on photosynthesis, Nfact is the nitrogen stress factor, Prft is the temperature stress factor, and PRFO3 is the photosynthetic reduction factor resulting from the interaction of ozone stress with water and atmospheric CO2 stress. To find the minimum value of the function.

[0034] The formula for calculating the leaf senescence acceleration factor caused by ozone stress is as follows: Among them, SLFO3 is the leaf senescence acceleration factor caused by ozone stress, SLFOZ1 and SLFOZ2 are set to 0.008 and 0.800 (tolerant type), 0.025 and 0.375 (moderately sensitive type) and 0.040 and 0 (sensitive type) respectively according to the wheat's tolerance to ozone, and OZONX is the target ozone concentration index.

[0035] The fourth sub-step is the field management document (.WHX).

[0036] In this embodiment, field management data serves as a crucial foundation for simulation accuracy assessment and parameter localization. This data primarily includes management information such as land preparation methods, sowing dates, sowing methods, sowing depth, sowing density, irrigation dates, irrigation amounts, fertilization dates, fertilizer types, and fertilizer amounts, as well as growth and yield records for winter wheat, including phenological indicators, biomass, grain weight, thousand-grain weight, and actual yield. Winter wheat monitoring data from agricultural meteorological stations largely meets these requirements. For instances where management data such as fertilizer types, amounts, and irrigation amounts are missing from certain stations in specific years, similar data from neighboring years or multiple years are used to supplement the data. It should be noted that if the winter wheat yield per unit area recorded by an agricultural meteorological station shows significant anomalies, then neighboring field observation records or average yield records from the county where the station is located are used to supplement the data.

[0037] The second step is to determine the optimal set of model parameters for the O3 modified-NWheat model.

[0038] Among them, the ozone stress function was used to improve the NWheat model, and the leaf senescence rate function was determined to obtain the O3 modified-NWheat model.

[0039] It should be noted that the NWheat model (wheat growth model) not only simulates the growth and development of wheat, water and nitrogen cycles, and the effects of various forcing factors on wheat growth and development on a diurnal scale, but also determines the impact of increased atmospheric CO2 on crop growth, radiation use efficiency, and transpiration efficiency based on CO2, radiation use efficiency, and temperature functions (CO2 and transpiration efficiency functions). Photosynthesis and leaf senescence depend on the duration and severity of abiotic stress; therefore, NWheat employs two different subroutines to calculate the effects of abiotic stress on reduced photosynthesis and accelerated leaf senescence.

[0040] In this embodiment, the ozone stress function incorporates the effects of ozone on wheat into the NWheat model by modifying the reduction in photosynthesis and the rate of leaf senescence. This allows the model to account for the influence of ozone in CO2, radiation use efficiency, temperature function, radiation use efficiency, and water stress function. The ozone stress function in the NWheat model mainly includes two reactions: first, ozone directly affects photosynthesis through stomatal flux; second, ozone accelerates leaf senescence and indirectly affects photosynthesis. These two reactions have different degrees of influence and are affected by other growth environments. For example, water stress or increased atmospheric CO2 concentration will lead to stomatal closure, in which case the direct impact of ozone on photosynthesis will be affected by the stomatal flux. Based on the existing abiotic stress methods in the NWheat model, the photosynthetic reduction effect of ozone stress is also incorporated into the model, and the interaction between ozone concentration and water stress and increased atmospheric CO2 concentration is considered simultaneously. The accelerating effect of ozone on leaf senescence is incorporated into a process similar to heat stress; ozone concentrations above a set threshold will lead to accelerated leaf senescence.

[0041] After modifying the NWheat model using the ozone stress function, the O3 modified-NWheat model was determined, and the leaf senescence rate function for each subregion can be obtained. The calculation formula is as follows: Among them, Leafsen is the daily total leaf senescence rate, Slan is the normal leaf senescence level based on the growth stage after emergence, Sfactor is the leaf senescence level caused by abiotic stress (mainly temperature and water), and SLFO3 is the leaf senescence acceleration factor caused by ozone stress.

[0042] Thus, this embodiment obtains an O3 modified-NWheat model affected by ozone.

[0043] Because the initial development environment of the model differs from the actual development environment, the model parameters should be localized according to the environment of the study area and the characteristics of the crop variety before being used in a different location. This involves obtaining an optimal set of model parameters, which is a prerequisite for model usability and simulation performance. The files that determine the model parameters include species parameter files (.SPE files), ecotype parameter files (.ECO files), and variety genetic parameter files (.CUL files). Species parameters depend on the crop itself, while ecotype parameters are mainly related to the crop's growth environment and environmental sensitivity. Generally, these are not adjusted for the crop's growth environment and environmental sensitivity; therefore, model parameter localization primarily involves processing the genetic parameters of the crop variety.

[0044] Currently, model parameter localization methods mainly fall into two categories. The first, the frequentist approach, estimates parameters based on data samples. By setting fixed parameters and creating an array, the estimated parameter values ​​are made closer to the true values. However, this method ignores the prior distribution information of the parameters, resulting in low accuracy of the localized parameters. The second approach, Bayesian estimation, considers both data samples and prior parameter information. This method takes into account the posterior probability distribution of the parameters, improving the accuracy of the localized parameters and gaining wider application. This embodiment employs the GLUE (Generalized Likelihood Uncertainty Estimation) method for crop genetic parameter localization based on Bayesian theory. First, it generates a large number of parameter values ​​through prior distribution to discretize the parameter space. Then, it calculates the likelihood value of each parameter based on field measurement data. Finally, it uses the Bayesian algorithm to calculate the posterior probability distribution of the parameters. The specific implementation steps include: The first sub-step involves obtaining the first preset number of variety genetic parameters in the variety genetic parameter file based on the input data required by the crop state estimation model.

[0045] In this embodiment, to verify the portability of localized parameters, the parameter data for historical periods are divided into training parameters and validation parameters. Data from 2015 to 2016 is used for parameter validation, and data from 2017 to 2018 is used for model performance validation. First, based on the input data required by the crop state estimation model, soil files, meteorological files, and field management files are created. A first preset number of variety genetic parameters is created in the variety genetic parameter file. The first preset number can be 25. Each set of variety genetic parameters includes 8 crop genetic parameters and 3 soil phenotypic parameters. The minimum number of GLUE parameter runs required for adjustment is 6000. To increase the parameter sample size, the number of GLUE runs is set to 20000 during parameter optimization.

[0046] It should be noted that this embodiment mainly calibrates and validates the model based on four key crop phenotypes: flowering date, maturity date, grain number, and grain weight. Eight varietal genetic parameters were selected for localization: VSEN, PPSEN, P1, P5, PHINT, GRNO, MXFIL, and STMMX. Among these, VSEN, PPSEN, P1, and PHINT are mainly related to flowering date; P5 is mainly related to maturity date; GRNO is mainly related to grain number; and MXFIL and STMMX are mainly related to grain weight. The physical meaning and parameter range of each varietal genetic parameter are shown in Table 1. Table 1 Since the uncertainty brought about by soil parameters is often considered to be greater than the uncertainty of the model itself, three soil phenotypic parameters were optimized based on the sensitivity analysis: soil runoff curve number (SLRO, 50~100), soil mineralization factor (SLNF, 0~1), and soil photosynthesis factor (SLPF, 0~1).

[0047] It is worth noting that regarding the parameter localization principle, during the parameter localization process, it is necessary to adhere to the principle of first checking parameters related to the critical growth period, and then checking parameters related to yield. During parameter localization, based on ensuring that the error between the model's simulated critical growth period and the observed value is controlled within 7 days, yield-related parameters are then screened, and finally, the set of parameters with the best yield simulation performance is selected as the model parameters.

[0048] The second sub-step involves using the GLUE method to optimize the genetic parameters of each group of varieties until the GLUE is run a fourth preset number of times, thereby obtaining the model output data corresponding to the optimized genetic parameters of each group of varieties.

[0049] In this embodiment, the fourth preset number can be 20,000. By performing 20,000 parameter simulations on each group of variety genetic parameters, the optimized model output data corresponding to each group of variety genetic parameters can be obtained. The model output data is the simulated crop state value. The implementation process of the GLUE method is existing technology and is not within the scope of protection of this invention; therefore, it will not be described in detail here.

[0050] The third step involves calculating the weights, mean, and variance of the genetic parameters for each variety based on the simulated and observed crop state values ​​corresponding to each group of genetic parameters. Then, based on these weights, mean, and variance, the optimal model parameter set for the crop state estimation model is determined. This process includes the following steps: Based on the simulated and observed crop state values ​​corresponding to the genetic parameters of each variety, the likelihood value between the simulated and observed crop state values ​​is calculated.

[0051] In this embodiment, to facilitate subsequent analysis of the simulation performance of each set of parameters, it is necessary to calculate the likelihood value between the simulated and observed values ​​of each set of parameters. The likelihood value between two sets of values ​​can be the mean of the corresponding distances between the two sets of values; that is, the mean of the distances is determined as the likelihood value between the two sets of values. Similarly, the correlation coefficient between the two sets of data can be calculated, and the final correlation coefficient can be used as the final likelihood value.

[0052] Based on the likelihood values ​​corresponding to the genetic parameters of each group of varieties, calculate the weight, mean, and variance of the genetic parameters for each group of varieties. The formula for calculating the weight of the genetic parameters for each group of varieties is as follows: in, Let be the weights of the genetic parameters of the i-th variety group. Let be the likelihood value corresponding to the genetic parameters of the i-th variety group, N be the number of the genetic parameter groups of the variety, and i be the sequence number of the genetic parameter group of the variety.

[0053] The formula for calculating the mean of genetic parameters for each group of varieties is: in, The mean values ​​of the genetic parameters for each group of varieties are given. The weights of the likelihood values ​​corresponding to the genetic parameters of the i-th variety group are: Let N be the genetic parameters of the i-th variety group, N be the number of the variety genetic parameter groups, and i be the sequence number of the variety genetic parameter group.

[0054] The formula for calculating the variance of the genetic parameters for each group of varieties is as follows: in, The variance of the genetic parameters for each group of varieties is given. The weights of the likelihood values ​​corresponding to the genetic parameters of the i-th variety group are: Let N be the genetic parameters of the i-th variety group, and N be the number of variety genetic parameter groups. This represents the mean of the genetic parameters for each group of varieties.

[0055] Based on the parameter optimization results, the top 10% of data were selected according to the likelihood value, weight, mean, and variance. The genetic parameter set of the varieties with good simulation effects on key growth period and yield was determined as the optimal parameter set. Based on this, the winter wheat production in 2017-2018 was simulated to verify the transferability of the selected parameter set. Finally, the optimal parameter set was selected and written into the variety parameter file.

[0056] S3. Based on the O3 modified-NWheat model corresponding to the optimal model parameter set and the basic data files of each sub-study area, determine the potential yield loss caused by ozone to each sub-study area.

[0057] The basic data files of each sub-study area are used as input data and input into the O3 modified-NWheat model corresponding to the optimal model parameter set to determine the crop yield of each sub-study area under the influence of ozone. The O3 modified-NWheat model is then used to simulate the crop yield when it is not affected by ozone. Based on the crop yield of each sub-study area under the influence of ozone and the crop yield when it is not affected by ozone, the difference between the crop yield when it is not affected by ozone and the crop yield when it is affected by ozone is determined as the potential yield loss caused by the corresponding sub-study area.

[0058] In this embodiment, based on localized model parameters, the O3 modified-NWheat model was used to simulate the potential changes in winter wheat yield under actual conditions where ozone is absent (or the ozone concentration is below the 25 ppb threshold). The changes in winter wheat yield in major winter wheat producing areas from 2015 to 2018 after removing the influence of ozone are shown below. Figure 5 As shown. By Figure 5It is known that the yield loss per unit area of ​​winter wheat caused by ozone in the study area ranges from 0 to 984 kg / ha. The largest yield reduction occurred in the HBAC region in 2018, accounting for 22.16% of the region's yield per unit area that year. In contrast, the HBAD and HBAJ regions, due to their consistently low ozone concentrations, generally experienced yield losses of less than 1% due to ozone. The potential yield loss of winter wheat is highly correlated with the distribution trend of ozone concentration; areas with high yield reductions are also often high-yielding winter wheat areas, while areas with low yield reductions are mainly distributed in coastal areas and the southern part of the study area where ozone concentrations are relatively low. The statistical results of the potential yield loss of winter wheat caused by ozone in 25 regions within the main winter wheat producing area are as follows: Figure 6 As shown, by Figure 6 It can be seen that the potential yield loss caused by ozone was the largest in 2017, followed by 2018 and 2017, with the smallest yield loss in 2015. Between 2015 and 2018, the maximum potential yield loss of winter wheat due to ozone in my country's main winter wheat producing areas reached 621, 841, 737, and 984 kg / ha, respectively, with average values ​​of 306.92, 346.52, 454.56, and 418.50 kg / ha, respectively, accounting for 4.96%, 5.73%, 7.58%, and 7.33% of the average winter wheat yield in the region that year. In the core winter wheat producing areas (yield per unit area > 6000 kg / ha), the average yield loss caused by ozone over the four years was 376.47, 423.00, 568.53 and 492.45 kg / ha, respectively, accounting for 5.52%, 6.20%, 8.61% and 7.68% of the average yield of the corresponding regions in those years. This also shows that ozone has a greater impact on high-yield winter wheat areas in my country.

[0059] Thus, this embodiment accurately assessed the potential impact of ozone on winter wheat yield. Specifically, based on the NWheat crop model incorporating the mechanism of ozone's influence on wheat growth, the potential yield loss of winter wheat caused by ozone in major winter wheat producing areas between 2015 and 2018 was assessed. The results of model training and validation show that the O3 modified-NWheat model has good simulation performance for the key growth stages of winter wheat, achieving high simulation accuracy even without extensive parameter optimization.

[0060] Since simulation accuracy is a direct criterion for judging the usability of a model, we select commonly used simulation accuracy criteria to examine the simulation performance of the NWheat model, including the coefficient of determination (R²). 2 The indicators include linear regression equation, root mean square error (RMSE), standardized root mean square error (NRMSE), and mean deviation error (MBE). The calculation formulas for each indicator are as follows: The formula for calculating the coefficient of determination is as follows: in, As the coefficient of determination, For the observed values, These are simulated values. The average of the observed values. is the average of the simulated values, and n is the total number of samples.

[0061] Linear regression coefficients.

[0062] The linear regression coefficients (slope a and intercept b) between observed and simulated values ​​directly reflect the degree of consistency between them. When a approaches 1 and b approaches 0, the simulated and measured values ​​are close to being consistent. The calculation process of the linear regression coefficients is existing technology and is not within the scope of protection of this invention; therefore, it will not be described in detail here.

[0063] The formulas for calculating the root mean square error and the standardized root mean square error are as follows: in, The root mean square error, For the observed values, These are simulated values. To standardize the root mean square error, The maximum value of the observed value. Let n be the minimum value of the observed values, and n be the total number of samples.

[0064] The formula for calculating the mean offset error is as follows: in, This is the mean offset error. For the observed values, These are simulated values, and n is the total number of samples.

[0065] It should be noted that among the above judgment indicators, R... 2 The goodness-of-fit index is an important statistical indicator; the closer its value is to 1, the better the fit between the observed and simulated values. RMSE and NRMSE reflect the magnitude of the absolute error between the simulated and observed values. An NRMSE value less than 0.1 is often considered to indicate excellent simulation accuracy, while values ​​between 0.1 and 0.2, 0.2 and 0.3, and >0.3 are considered to indicate good, moderate, and poor simulation accuracy, respectively. MBE reflects the degree to which the model overestimates or underestimates the simulated object; a positive MBE indicates overestimation, and a negative MBE indicates underestimation, with a larger absolute value indicating greater overestimation or underestimation.

[0066] In most sub-regions of the main winter wheat producing areas, winter wheat yields can reach over 6000 kg / ha. High-yield areas are concentrated in Hebei, eastern Henan, western Shandong, and northern Anhui and Jiangsu. Lower-yield areas are mostly in southwestern Henan and southern Anhui and Jiangsu. These areas are also transitional zones between wheat and rice cultivation in my country. The spatial distribution of winter wheat yields recorded by 25 agricultural meteorological stations in the main winter wheat producing areas is shown below. Figure 7 As shown.

[0067] The flowering period is one of the most critical growth stages affecting winter wheat yield. The training results during the flowering period are shown in the figure below. Figure 8 The verification results during the flowering period are shown in the figure below. Figure 9 As shown, this indicates that NWheat has a good simulation effect on the flowering period, and the training and validation R values ​​are [data missing]. 2 The values ​​were 0.98 and 0.97, respectively, and the NRMSE values ​​were 0.04 and 0.05, indicating that the model achieved excellent simulation accuracy (NRMSE < 0.1). Furthermore, the RMSE for training and validation were 3.18 days and 3.99 days, respectively, and the MBE was -1.88 days and -1.71 days, respectively, indicating that the NWheat model simulated the flowering period earlier than the actual recorded flowering period, but the advance was less than 2 days in both cases. The spatial distribution of the error between the observed and simulated flowering periods of winter wheat in the study area is shown in the figure. Figure 10 Within the study area, the errors of both simulations and observations were within 7 days, and the error of the simulated flowering period was less than 3 days in most areas. Except for 2018, when the simulated flowering period was later than the observed value at most stations, the simulated flowering period in other years was mostly earlier than the actual flowering period.

[0068] The training results during the mature stage are shown in the figure below. Figure 11 The verification results at the maturity stage are shown in the figure below. Figure 12 As shown, R for model training and validation 2 The RMSE values ​​were 0.97 and 0.97 respectively, and the NRMSE values ​​were 0.05 and 0.04 respectively, indicating that the model achieved excellent simulation accuracy for the maturity period, similar to that for the flowering period. Unlike the simulation of the flowering period, the simulated maturity period was slightly delayed compared to the observed actual maturity period, with RMSE values ​​of 3.18 days and 3.46 days for training and 0.84 days and 1.40 days for validation. The spatial distribution of the error between the observed and simulated values ​​of winter wheat maturity in the study area is shown below. Figure 13 This indicates that, similar to the simulation results for the flowering period, the simulation error for the maturity period is also within 7 days, and except for 2018, the simulation error for the maturity period in most regions is less than 3 days. Overall, except for a few regions, the simulated maturity date for winter wheat is mostly later than the actual maturity date.

[0069] The training results of the model for simulating winter wheat yield are shown in the figure below. Figure 14 The model's validation results for winter wheat yield simulation are shown in the figure below. Figure 15 As shown, the training group's simulation of winter wheat yield was significantly better than the validation group, with R0... 2 The RMSE values ​​were 0.93 and 0.88, respectively, and the RMSE values ​​for the training and validation years were 272.35 and 402.93 kg / ha, respectively. This further illustrates that the training effect of the model is significantly better than the validation effect. Furthermore, the simulated winter wheat yield in the training group was slightly higher than the actual yield (MBE of 36.36 kg / ha), while the validation group showed a significant underestimation, with an overall underestimation of -163.46 kg / ha. Although the simulation effect of the validation group was slightly worse than that of the training group, the NRMSE values ​​were 0.08 and 0.09, respectively, which also indicates that the model achieved excellent simulation results for winter wheat yield. The spatial distribution of the error between the observed and simulated values ​​of winter wheat yield in the study area is shown in the figure. Figure 16 As shown, this indicates that the error between the simulated yield and the actual yield is less than 300 kg / ha in approximately 80% of the regions. Regions with yield errors greater than 600 kg / ha are mainly distributed in southern Shanxi, northern Henan, Hebei, and northern Jiangsu. Only a few regions have simulated yield errors reaching 1000 kg / ha, such as the HBAW region in 2017 and the HBAG region in 2018, with simulation errors of 1359 and 1231 kg / ha respectively, both of which are underestimates.

[0070] In summary, the O3 modified-NWheat model achieved excellent simulation performance for the key growth stages and yield of winter wheat in North China. The simulated growth period from flowering to maturity was slightly longer than the actual growth period, and the simulated yield was slightly lower than the actual yield. Therefore, the O3 modified-NWheat model can be used to assess the impacts of ozone and climate change on the growth process and yield of winter wheat in the study area. The localization results of winter wheat varieties and model parameters in the study area are shown in Table 2. Table 2 This invention provides a method for assessing the impact of surface ozone on winter wheat yield based on multi-factor stress. This method is based on the ozone-modified O3-NWheat model, which takes into account the combined stress effects of ozone and CO2, while considering different regional varieties, meteorological factors, soil factors, and field management factors. Taking the winter wheat growth and yield observation data of agricultural meteorological stations in the main winter wheat producing areas of my country as an example, the method quantitatively assesses the impact of ozone on winter wheat yield in my country after verifying the applicability of the model in major agricultural areas, and accurately assesses the degree of loss of potential crop yield caused by ozone.

[0071] The above-described embodiments are only used to illustrate the technical solutions of the present invention, and are not intended to limit it. Although the present invention has been described in detail with reference to the foregoing embodiments, those skilled in the art should understand that modifications can still be made to the technical solutions described in the foregoing embodiments, or equivalent substitutions can be made to some of the technical features. Such modifications or substitutions do not cause the essence of the corresponding technical solutions to deviate from the scope of the technical solutions of the embodiments of the present invention, and should all be included within the protection scope of the present invention.

Claims

1. A method for assessing the impact of surface ozone on winter wheat yield based on multi-factor stress, characterized in that, Includes the following steps: Based on the location of crop meteorological stations within the study area to be evaluated, the study area to be evaluated is divided into sub-study areas within the study area to be evaluated. Obtain the basic data files for each sub-study area and determine the optimal model parameter set corresponding to the O3 modified-NWheat model; Based on the O3 modified-NWheat model corresponding to the optimal model parameter set and the basic data files of each sub-study area, the potential yield loss caused by ozone in each sub-study area is determined. The basic data files include variety files, soil files, meteorological files, and field management files. The variety files are winter wheat varieties from each sub-study area during the study period. The soil files include station latitude and longitude, soil type, gravel content, soil color, surface slope, organic carbon content, pH value, soil moisture content, saturated moisture content, soil saturated hydraulic conductivity, field water holding capacity, and wilting point water holding capacity. The meteorological files include daily maximum temperature, daily minimum temperature, precipitation, solar radiation, wind speed, relative humidity, CO2 concentration, evapotranspiration, and target ozone index. The field management files include management data and growth and yield record data. The management data includes land preparation method, sowing date, sowing method, sowing depth, sowing density, irrigation date, irrigation amount, fertilization date, fertilizer type, and fertilizer amount. The growth and yield record data includes winter wheat phenological indicators, biomass, grain weight, thousand-grain weight, and actual yield. Determine the optimal set of model parameters for the O3 modified-NWheat model, including: Based on the input data required for the O3 modified-NWheat model, create a first preset number of variety genetic parameters in the variety genetic parameter file; wherein, a set of variety genetic parameters includes a second preset number of crop genetic parameters and a third preset number of soil phenotypic parameters; The GLUE method was used to optimize the genetic parameters of each group of varieties until the GLUE was run a fourth preset number of times, and the model output data corresponding to the optimized genetic parameters of each group of varieties was obtained; wherein, the model output data is the simulated value of crop status in each sub-study area; Based on the simulated and observed crop state values ​​corresponding to the genetic parameters of each variety, the weight, mean, and variance of the genetic parameters of each variety are calculated. Then, based on the weight, mean, and variance, the optimal model parameter set corresponding to the O3 modified-NWheat model is determined. Based on the optimal model parameter set corresponding to the O3 modified-NWheat model and the basic data files of each sub-study area, the potential yield loss caused by ozone in each sub-study area is determined, including: The basic data files of each sub-study area are used as input data and input into the O3modified-NWheat model corresponding to the optimal model parameter set to determine the crop yield of each sub-study area when affected by ozone. The O3modified-NWheat model is then used to simulate the crop yield when not affected by ozone. Based on the crop yield of each sub-study area when affected by ozone and when not affected by ozone, the difference between the crop yield when not affected by ozone and the crop yield when affected by ozone is determined as the potential yield loss caused by ozone in the corresponding sub-study area.

2. The method for assessing the impact of surface ozone on winter wheat yield based on multi-factor stress as described in claim 1, characterized in that, The process of calculating the weights, mean, and variance of the genetic parameters for each group of varieties includes: Based on the simulated and observed crop state values ​​corresponding to the genetic parameters of each variety, calculate the likelihood value between the simulated and observed crop state values. Based on the likelihood values ​​corresponding to the genetic parameters of each group of varieties, calculate the weight, mean, and variance of the genetic parameters of each group of varieties.

3. The method for assessing the impact of surface ozone on winter wheat yield based on multi-factor stress as described in claim 2, characterized in that, The formula for calculating the weight of the genetic parameters for each group of varieties is as follows: in, Let be the weights of the genetic parameters of the i-th variety group. Let N be the likelihood value corresponding to the genetic parameters of the i-th variety, N be the number of the genetic parameter groups of the variety, and i be the sequence number of the genetic parameter group of the variety. The formula for calculating the mean of genetic parameters for each group of varieties is: in, The mean values ​​of the genetic parameters for each group of varieties are given. The weights of the likelihood values ​​corresponding to the genetic parameters of the i-th variety group are: Let N be the genetic parameters of the i-th variety group, N be the number of the variety genetic parameter groups, and i be the sequence number of the variety genetic parameter group. The formula for calculating the variance of genetic parameters for each group of varieties is: in, The variance of the genetic parameters for each group of varieties is given. The weights of the likelihood values ​​corresponding to the genetic parameters of the i-th variety group are: Let N be the genetic parameters of the i-th variety group, and N be the number of variety genetic parameter groups. This represents the mean of the genetic parameters for each group of varieties.

4. The method for assessing the impact of surface ozone on winter wheat yield based on multi-factor stress as described in claim 1, characterized in that, Based on the location of crop meteorological stations within the study area to be evaluated, the study area is divided into sub-study areas, including: Using the Thiessen polygon method, the study area to be evaluated is divided into sub-study areas based on the location of crop meteorological stations within the study area to be evaluated. Among them, the edge distance of the sub-study area is the shortest straight-line distance to the corresponding crop meteorological station.

Images