Field-level rice yield remote sensing estimation method based on high-resolution total primary productivity and dynamic harvest index

Abstract

The invention provides a field-level rice yield remote sensing estimation method based on high-resolution total primary productivity and a dynamic harvest index, and the method comprises the following steps: firstly, obtaining a rice growth season time sequence Sentinel-2 image based on a remote sensing cloud platform, calculating a leaf area index, and obtaining meteorological data EAR5 products and CO2 data; then, a self-service sampling method and a multi-model process sensitivity analysis method are used for analyzing the influence rule of the rice growth season environment factors on the GPP, an improved double-leaf light energy utilization rate model is constructed, and high-resolution rice GPP high-precision estimation is achieved. And finally, forming an HIDYM method in combination with the dynamic harvest index to realize the yield estimation of the field-level rice in Jiangsu province for many years. The method can be used for quickly and accurately estimating the yield of the large-range field-level rice.

Description

Technical Field

[0001] The invention belongs to the field of agricultural remote sensing monitoring, and is a field-level rice yield remote sensing estimation method based on high-resolution gross primary productivity and dynamic harvest index. Background Art

[0002] In recent decades, domestic and foreign scholars have developed a series of satellite remote sensing estimation methods for crop yield, mainly including empirical statistics / data-driven methods, semi-empirical and semi-mechanistic methods based on gross primary productivity (GPP), and crop growth model data assimilation methods. These methods generally have problems such as low accuracy, low efficiency, or poor universality. Among them, the semi-empirical and semi-mechanistic method based on GPP has great potential in large-scale crop yield remote sensing prediction. It converts aboveground biomass into grain yield through harvest index based on the consideration of photosynthesis, dry matter accumulation and distribution physiological mechanisms. They use a large number of satellite remote sensing images, not only showing the ability to capture the spatiotemporal pattern of crop yield, but also being able to be applied in new regions and years without relying on ground-based measured yield data for recalibration. Many studies have demonstrated the ability of semi-empirical and semi-mechanistic models to draw spatial distribution maps of crop yield using county-level statistical yield data. However, it is still unclear whether semi-empirical and semi-mechanistic models can accurately estimate the yield of individual plots that are above or below the county average. Therefore, it is necessary to evaluate the performance of semi-empirical and semi-mechanistic models over a large area (province / state) using field-measured yield data, and to study the spatiotemporal variability of crop yield using the increasingly abundant Sentinel-2 images. These high-resolution spatial distribution maps of yield will provide us with new insights into the impact of climate change, crop management practices, and pests and diseases on crop production.

[0003] As the main variable in the semi-empirical and semi-mechanistic yield estimation model, GPP is usually estimated by sampling statistical models, process models or light energy utilization models. In particular, the light energy utilization model driven by remote sensing data is widely used in regional GPP estimation due to its simplicity and strong physiological mechanism. However, there are still problems such as large uncertainty in the estimation of GPP in farmland ecosystems and low spatial resolution. In addition, most of the current crop yield estimation studies based on the semi-empirical and semi-mechanistic model of GPP use a fixed value of the harvest index of a given crop across years or regions. However, as an important indicator sensitive to crop categories, crop varieties, field management measures and environmental conditions, the harvest index often shows significant differences in different years and regions, and has obvious characteristics of temporal and spatial variation. Previous studies have shown that when GPP is accurately estimated (about 90%), the estimated yield deviates greatly from the measured value. The main reason is that these models assume that the harvest index is a constant and do not consider the impact of environmental stress on the distribution of crop dry matter among different organs. Therefore, there is an urgent need to build an automated method based on high-resolution GPP and dynamic harvest index for efficient estimation and mapping of rice grain yield in large areas of farmland without the need for model training through a large number of field measured yield samples. Summary of the invention

[0004] The technical problem solved by the present invention is to provide a field-level rice yield remote sensing estimation method based on high-resolution gross primary productivity and dynamic harvest index. First, the leaf area index LAI is obtained by inverting the hybrid model of the PROSAIL model and the neural network algorithm using the rice growing season time series Sentinel-2 image, determining the start and end dates of the rice growing season, and obtaining the daily average temperature and other variables of the rice growing season meteorological data product ERA5-Land; then, the original two-leaf light energy utilization model is transformed into 5 variants by the self-service sampling method, and the GPP accuracy is verified based on the flux tower; then, the multi-model process sensitivity analysis method xSSA is used to analyze the influence of environmental factors on GPP during the rice growing season, and an improved two-leaf light energy utilization model mTL-LUE is constructed, thereby realizing high-precision estimation of high spatial resolution GPP and growing season GPP accumulation. Finally, by coupling the dynamic harvest index estimated by the PDR-TCT method, a rice yield estimation method HIDYM with high-altitude resolution GPP and dynamic harvest index is formed, and a multi-year large-scale 10m resolution rice yield spatial distribution product is generated. At the same time, the performance of HIDYM was systematically evaluated using years of field-measured yield data and compared with the widely used satellite remote sensing GPP estimation product GPP MOD17The method is compared with the fixed harvest index strategy. This method has a certain theoretical basis, fewer parameters and a simple model structure. It is easy to migrate to other years and regions, and has a strong advantage in estimating crop yields in areas with complex vegetation types, large areas and long periods of time.

[0005] The technical solution to achieve the purpose of the present invention is:

[0006] A method for estimating rice yield at field level by remote sensing based on high-resolution gross primary productivity and dynamic harvest index comprises the following steps:

[0007] Step 1: Obtain Sentinel-2 satellite images of the rice growing season in the study area, preprocess and cloud mask the images on the remote sensing cloud platform, and calculate the leaf area index LAI based on the PROSAIL model and the neural network hybrid model. The Sentinel-2 image calculation uses 8 bands of information. Finally, the rice growing season time series LAI data and the inflection point method are used to calculate the start and end dates of the rice growing season.

[0008] Step 2: Obtain the daily average temperature T during the rice growing season based on the meteorological data ERA5-Land product of the study area a , relative humidity RH and dew point temperature, calculate the saturated water vapor pressure difference VPD; at the same time, use the downlink shortwave radiation to calculate the photosynthetically active radiation PAR, in addition, the atmospheric CO2 data is obtained synchronously at the Earth System Research Laboratory of the National Oceanic and Atmospheric Administration in the United States;

[0009] Step 3: The original two-leaf light energy utilization model was transformed into five forms using the self-service sampling method, taking into account the atmospheric CO2 fertilization effect and environmental factors. The original model and five variants were evaluated for accuracy based on the gross primary productivity (GPP) measured by the flux tower to determine the optimal rice GPP estimation model.

[0010] Step 4: To clarify the environmental factors VPD, T a In order to investigate the physical mechanism of the nonlinear effect of CO2 on GPP, this study used the multi-model process sensitivity analysis method xSSA to analyze the GPP AND 、GPP OR1 、GPP OR2 、GPP OR3 and GPP OR4 Conduct sensitivity analysis of the model and model parameters to assess the impact of model inputs on the uncertainty of GPP estimates and provide sensitivity measures to process options and model processes;

[0011] Step 5: Based on the optimal rice GPP estimation model determined by the bootstrap sampling method and multi-model process sensitivity analysis method in steps 3 and 4, an improved two-leaf light energy utilization model was constructed to estimate the GPP of rice in all fields in Jiangsu Province from 2019 to 2022.

[0012] In a preferred embodiment, the present invention further comprises step 6: coupling the high-resolution GPP accumulation amount and the dynamic harvest index during the rice growing season in step 5 to form a rice yield estimation method HIDYM driven by high spatial resolution GPP and dynamic harvest index, so as to achieve rice yield estimation for all fields in Jiangsu Province from 2019 to 2022. At the same time, the performance of HIDYM was systematically evaluated using field measured yield data for many years, and compared with the widely used satellite remote sensing GPP estimation product GPP in terms of rice and winter wheat yield prediction in Jiangsu Province. MOD17 It is compared with the fixed harvest index strategy.

[0013] Furthermore, the field-level rice yield remote sensing estimation method based on high-resolution gross primary productivity and dynamic harvest index proposed in the present invention comprises the following specific steps in step 1:

[0014] Step 1-1: Obtain Sentinel-2 images of the rice growing season in the study area on the remote sensing cloud computing platform and perform cloud mask processing;

[0015] Step 1-2: Use the PROSAIL model and the neural network hybrid model to invert the 8-band reflectance information of the Sentinel-2 image to obtain the rice leaf area index LAI;

[0016] Step 1-3: Use the rice growing season time series LAI and the inflection point method to calculate the start and end dates of the rice growing season.

[0017] Furthermore, the field-level rice yield remote sensing estimation method based on high-resolution gross primary productivity and dynamic harvest index proposed in the present invention, step 2 specifically comprises the following steps:

[0018] Step 2-1: Obtain the daily average temperature T of the ERA5-Land meteorological data product during the rice growing season in the study area a , relative humidity RH, dew point temperature and downward shortwave radiation;

[0019] Step 2-2: At the same time, obtain atmospheric CO2 data from the Earth System Research Laboratory of the National Oceanic and Atmospheric Administration of the United States;

[0020] Step 2-3: Calculate the saturated vapor pressure difference VPD using the average air temperature, relative humidity and dew point temperature, and calculate the photosynthetically active radiation PAR using the downward shortwave radiation;

[0021] Step 2-4: Use the coordinates of the Jurong Ecological Experiment Station in Jiangsu Province in the previous n years to extract the above meteorological indicators and the remote sensing indicators in step 1.

[0022] Furthermore, the field-level rice yield remote sensing estimation method based on high-resolution gross primary productivity and dynamic harvest index proposed in the present invention, the specific steps in step 3 include:

[0023] Step 3-1: Use the bootstrap sampling method to consider the atmospheric CO2 fertilization effect and environmental factors. Use logic to transform the original two-leaf light energy utilization model into 5 forms of GPP AND 、GPP OR1 、GPP OR2 、GPP OR3 and GPP OR4 Model;

[0024] Step 3-2: Compare the accuracy of the original two-leaf solar energy utilization model and the five variants with the GPP measured by the flux tower in different years to identify the optimal model for GPP estimation;

[0025] Step 3-3: Compare the daily changes of the original two-leaf light energy utilization efficiency model and the five variants with the GPP measured by the flux tower during the rice growing season in different years;

[0026] Step 3-4: Apply the coefficient of determination R 2 , root mean square error RMSE, and relative root mean square error RRMSE to evaluate and verify the results.

[0027] Furthermore, the field-level rice yield remote sensing estimation method based on high-resolution gross primary productivity and dynamic harvest index proposed in the present invention, the specific steps in step 4 include:

[0028] Step 4-1: Use multi-model process sensitivity analysis to analyze GPP AND 、GPP OR1 、GPP OR2 、GPP OR3 and GPP OR4 Sensitivity analysis was conducted using five models;

[0029] Step 4-2: Use multi-model process sensitivity analysis method to analyze model parameters (APAR sh ,APAR su ,g(CO2),g(T a ) and g(VPD)) were used to conduct global sensitivity analysis.

[0030] Furthermore, the field-level rice yield remote sensing estimation method based on high-resolution gross primary productivity and dynamic harvest index proposed in the present invention has the following specific steps in step 5:

[0031] Step 5-1: Determine the optimal model for GPP estimation and the most sensitive model structure and variables based on steps 3 and 4;

[0032] Step 5-2: Determine the optimal GPP estimation model as the improved two-leaf light energy utilization efficiency model mTL-LUE;

[0033] Step 5-3: Use the mTL-LUE model to estimate the cumulative GPP of rice growing season in Jiangsu Province in multiple years and draw the spatial distribution map of Jiangsu Province each year.

[0034] Furthermore, the field-level rice yield remote sensing estimation method based on high-resolution gross primary productivity and dynamic harvest index proposed in the present invention has the following specific steps in step 6:

[0035] Step 6-1: The high spatial resolution GPP accumulation during the rice growing season obtained in step 5 is converted into rice aboveground biomass through carbon utilization efficiency, carbon conversion to biomass coefficient and aboveground biomass allocation coefficient;

[0036] Step 6-2: The converted aboveground biomass is coupled with the PDR-TCT method to estimate the spatiotemporal dynamic index, forming a crop yield estimation method HIDYM driven by high spatial resolution GPP and dynamic harvest index;

[0037] Step 6-3: Use multiple years of field-measured yield data to systematically evaluate the performance of HIDYM and compare it with the widely used satellite remote sensing GPP estimation product GPP for rice yield estimation in Jiangsu Province. MOD17 Compare with the fixed harvest index strategy;

[0038] Step 6-4: Use the HIDYM method to estimate rice yield at the field level in Jiangsu Province for multiple years and generate a multi-year, large-scale, 10-meter-resolution rice yield spatial distribution product.

[0039] Step 6-5: Apply the coefficient of determination R 2 , root mean square error RMSE, and relative root mean square error RRMSE to evaluate and verify the results.

[0040] Preferably, the dynamic harvest index HI d Obtained by the following steps:

[0041] Step 601: Based on the multi-temporal Sentinel-2 satellite images within the three key time points of the beginning, peak and end of the rice growing season in the study area, the vegetation index and tassel cap transformation of all Sentinel-2 images within the three key time points are calculated; specifically, the following steps are included:

[0042] Step 601-1: Obtain information on rice planting distribution and rice growth process in the study area;

[0043] Step 601-2: Preprocessing and cloud masking of multi-temporal Sentinel-2 images to obtain cloud-free images;

[0044] Step 601-3: Calculate the Sentinel-2 vegetation index and tasseled cap transformation within the range of the start, peak and end dates of the rice growing season using the rice growth process information in the study area;

[0045] Step 602: respectively obtaining the vegetation index and tassel cap transformation components of the images at the time points corresponding to the start, peak and end dates of the rice growing season; specifically including:

[0046] Step 602-1: Extract the start, peak and end dates using the Sentinel-2 vegetation index (EVI) and inflection point method for the study area;

[0047] Step 602-2: Obtain the Sentinel-2 image vegetation index EVI and LSWI and the three components of tasseled cap transformation, namely brightness, humidity and greenness, for the start, peak and end dates of rice production;

[0048] Step 603: According to the definition of harvest index and the yield formation process, a phenological difference ratio method PDR is proposed based on Sentinel-2 images of the start, peak and end dates of rice, and then the vegetation index and tasseled cap transformation components are coupled to form a new remote sensing index. A dynamic harvest index estimation model is constructed in combination with the measured harvest index. Finally, the dynamic harvest index estimation model constructed by the vegetation index and tasseled cap transformation components is compared in accuracy by using the exhaustive method and feature screening method to determine the optimal dynamic harvest index estimation model; specifically, the following steps are performed:

[0049] Step 603-1: Obtain the phenological difference ratio PDR by the following formula:

[0050]

[0051] Among them, D SOS , D POS and D EOS vegetation index or tasseled cap transform representing the start, peak, and end dates of the rice growing season, respectively;

[0052] Step 603-2: coupling the phenological difference ratio PDR with the vegetation index VI and the tasseled cap transform TCT to form a remote sensing index with physiological significance, and then combining it with the measured harvest index to construct a dynamic harvest index estimation model;

[0053] Step 603-3: using exhaustive method and feature screening method to compare the accuracy of different combinations of vegetation index and tasseled cap transformation components, and determine the optimal dynamic harvest index estimation model;

[0054] Step 603-4: Apply the coefficient of determination R 2 , root mean square error RMSE, relative root mean square error RRMSE evaluation and verification results;

[0055] Step 604: For the land parcel to be estimated, obtain multi-temporal Sentinel-2 satellite images and input the optimal harvest index estimation model to obtain a dynamic harvest index.

[0056] Specific:

[0057] Vegetation index VI includes enhanced vegetation index EVI and land surface water index LSWI, and the calculation formula is as follows;

[0058]

[0059] where ρ NIR , ρ RED , ρ BLUE and ρ SWIR Represent the reflectance values ​​of Sentinel-2 band 8, band 4, band 2 and band 11 respectively; band 8 is 785-900nm, band 4 is 650-680nm, band 2 is 458-523nm, and band 11 is 1565-1655nm;

[0060] The three components of the Tasseled Cap Transform (TCT) are calculated using the coefficients of the 13 bands of the Sentinel-2 image. The specific calculation formula is as follows:

[0061]

[0062] Where m = 1, 2, 3, TCT1 represents the brightness component of the tasseled cap transform TCT, TCT2 represents the green component of the tasseled cap transform TCT, and TCT3 represents the humidity component of the tasseled cap transform TCT; B i Represents the reflectance value of the i-th band of the Sentinel-2 image; a i Represents the coefficient of the i-th band.

[0063] Preferably, in step 3-2, the remote sensing indicators PDR-TCT and PDR-VI with physiological significance have formulas as follows:

[0064]

[0065] Among them, TCT SOS 、TCT POS and TCT EOS TCT represents the start date, peak date and end date of the rice growing season m , m = 1, 2, 3, representing the brightness component, green component and humidity component of the tasseled cap transform TCT respectively;

[0066]

[0067] Among them, VI SOS VI POS and VI EOS The vegetation indices VI, which represent the start date, peak date, and end date of the rice growing season, include the Enhanced Vegetation Index EVI and the Land Surface Water Index LSWI.

[0068] Compared with the prior art, the present invention adopts the above technical solution and has the following technical effects:

[0069] 1. The field-level rice yield remote sensing estimation method based on high-resolution gross primary productivity and dynamic harvest index of the present invention is based on crop photosynthesis, dry matter accumulation and distribution physiological mechanisms, monitors the yield formation process, and has a strong agronomic mechanism;

[0070] 2. The field-level rice yield remote sensing estimation method based on high-resolution gross primary productivity and dynamic harvest index of the present invention has fewer parameters and a simple model structure, and is relatively easy to migrate to other years and regions. It has a strong advantage in estimating crop yields in areas with complex vegetation types, large areas, and long periods of time;

[0071] 3. The field-level rice yield remote sensing estimation method based on high-resolution gross primary productivity and dynamic harvest index of the present invention can achieve accurate and efficient estimation of crop yields over a large area and effectively reveal the spatiotemporal variability of rice yields in all fields across the province. BRIEF DESCRIPTION OF THE DRAWINGS

[0072] Figure 1 It is the technical roadmap of the field-level rice yield remote sensing estimation method based on high-resolution gross primary productivity and dynamic harvest index of the present invention.

[0073] Figure 2 The GPP (GPP) estimated by the TL-LUE model based on the field-level rice yield remote sensing estimation method based on high-resolution gross primary productivity and dynamic harvest index of the present invention is org ) and the GPP estimated based on the five variants of the mTL-LUE model and the GPP measured at the flux station (GPP EC ) on a daily scale.

[0074] Figure 3 The GPP estimated by different models of the field-level rice yield estimation method based on high-resolution gross primary productivity and dynamic harvest index of the present invention and the GPP measured at the flux station (GPP EC ) seasonal changes.

[0075] Figure 4Sobol' sensitivity analysis (A) of the mTL-LUE model of the field-level rice yield estimation method based on high-resolution total primary productivity and dynamic harvest index of the present invention. The five parameters of the model (APAR sh ,APARsu,g(CO2),g(T a ), g(VPD)), (B) Sobol' sensitivity index of five process options and (C) two processes (P1 and P2). P1 is the process; P2 is VPD, T a The Sobol' first-order sensitivity and total sensitivity are represented by dark and light bars, respectively.

[0076] Figure 5 This is the spatial distribution map of rice AGPP in Jiangsu Province from 2019 to 2022 based on the field-level rice yield remote sensing estimation method based on high-resolution gross primary productivity and dynamic harvest index, where (A) to (D) represent the years 2019-2022, respectively.

[0077] Figure 6 Is to use R 2 , RMSE and RRMSE indicators were used to evaluate the accuracy of rice yield estimation at the field level by coupling high-resolution GPP with dynamic harvest index or fixed harvest index. (A) and (B) represent the dynamic harvest index strategy and the fixed harvest index strategy, respectively.

[0078] Figure 7 It is a field-level rice yield remote sensing estimation method based on high-resolution gross primary productivity and dynamic harvest index for 2019-2022 based on GPP. OR3 and GPP MOD17 Comparison of model rice yield estimation accuracy.

[0079] Figure 8 It is a 10m spatial resolution rice yield spatial distribution product based on the field-level rice yield remote sensing estimation method using high-resolution gross primary productivity and dynamic harvest index, where (A) (B) (C) and (D) represent the years 2019-2022, respectively. DETAILED DESCRIPTION

[0080] The embodiments of the present invention are described in detail below, examples of which are shown in the accompanying drawings, wherein the same or similar reference numerals throughout represent the same or similar elements or elements having the same or similar functions. The embodiments described below with reference to the accompanying drawings are exemplary and are only used to explain the present invention, and cannot be interpreted as limiting the present invention.

[0081] The present invention is implemented based on time series Sentinel-2 satellite images and meteorological data. The input data set used to drive the model is shown in Table 1.

[0082] Table 1 Input dataset used to drive the mTL-LUE model

[0083]

[0084] This embodiment is carried out in Jiangsu Province, a major rice-growing province in China. The rice planting area in Jiangsu Province exceeds 33 million mu, and the output exceeds 19.8 million tons, ranking fourth in the country. It is mainly distributed in Yancheng City, Huai'an City, Nantong City, Taizhou City and other regions. In addition, the measured data of field-level rice yield in Jiangsu Province from 2019 to 2022 was used to conduct a systematic evaluation and verification of the HIDYM method, which can effectively verify the advancement, universality and robustness of the present invention.

[0085] A method for estimating rice yield at field level using remote sensing based on high-resolution gross primary productivity and dynamic harvest index. Figure 1 As shown in the figure, it mainly includes time series Sentinel-2 preprocessing and LAI inversion, ERA5-Land meteorological products and atmospheric CO2 data during the rice growing season, multi-year field-scale rice yield collection, self-service sampling method and multi-process model sensitivity analysis, construction of an improved two-leaf light energy utilization model, coupling dynamic harvest index to form HIDYM method, field-scale rice yield prediction and evaluation, and GPP MOD17 Compared with the fixed harvest index strategy, the rice yield mapping for different scales and years includes the following steps:

[0086] Step 1: Obtain the time series Sentinel-2 satellite images of the rice growing season in Jiangsu Province and preprocess the Sentinel-2 images on the remote sensing cloud platform. The preprocessing steps of Sentinel-2 include cloud removal, linear time interpolation and LAI inversion. The Sentinel-2 cloud probability dataset can characterize the probability that each pixel in the Sentinel-2 image is contaminated by clouds. The present invention first removes clouds from the original Sentinel-2 surface reflectance image based on this dataset to obtain a cloud-free Sentinel-2 image; then, in order to solve the data missing caused by cloud pollution, the linear time interpolation method is used on the GEE platform to fill in the missing values. The present invention calculates the rice LAI based on the reflectivity of the original 8 bands of Sentinel-2 and the PROSIAL model coupled neural network hybrid inversion algorithm. The accuracy of the inverted LAI is evaluated by comparing it with the measured rice LAI in Suining County. It is found that the root mean square error of the LAI inversion during the tillering, heading and filling stages is between 0.91 and 1.30 m 2 / m 2Subsequently, the 5-day LAI dataset was smoothed using the Savitzky-Golay filter and linear interpolation was performed to generate daily LAI products. Finally, the rice growing season time series LAI data and the inflection point method were used to obtain the start and end dates of the rice growing season.

[0087] Step 2: Calculate the meteorological indicators of the whole rice growth period according to the meteorological data of the study area and the start and end dates of the rice growing season, including:

[0088] Step 2-1: Use the start and end dates of the rice growing season calculated in step 1 to calculate the daily average temperature T during the entire rice growth period. a , relative humidity RH, daily dew point temperature, daily downward shortwave radiation;

[0089] Step 2-2: Calculate the saturated vapor pressure difference VPD using the daily average temperature, relative humidity and dew point temperature;

[0090] Step 2-3: Calculate daily photosynthetically active radiation PAR using daily downward shortwave radiation;

[0091] Steps 2-4: Synchronously acquire atmospheric CO2 data at the National Oceanic and Atmospheric Administration Earth System Research Laboratory.

[0092] Step 3: Use the bootstrap sampling method to obtain five variants of the original two-leaf light energy utilization model and determine the optimal GPP estimation model based on the flux tower data, including:

[0093] Step 3-1: Use the remote sensing and meteorological data obtained in steps 1 and 2 to drive the original two-leaf light energy utilization model;

[0094] Step 3-2: Considering the calculation rules of atmospheric CO2 fertilization effect and environmental limiting factors, the original two-leaf light energy utilization model was constructed into GPP using the self-service sampling method. AND 、GPP OR1 、GPP OR2 、GPP OR3 and GPP OR4 Five models,;

[0095] Step 3-3: Compare the accuracy of the original two-leaf solar energy utilization model and the five variants with the GPP measured by the flux tower in different years to identify the optimal model for GPP estimation ( Figure 2 ); As can be seen from the figure, compared with the GPP estimated based on the TL-LUE model (GPP org : R 2 =0.68,RMSE=2.89g C m -2 ), in the 5 variants of the mTL-LUE model, except for GPP OR4 (R 2=0.63,RMSE=3.15 g C m -2 ) In addition to GPP EC Showing a stronger correlation (GPP AND : R 2 =0.77, RMSE = 2.48 gC m -2 ;GPP OR1 : R 2 =0.72, RMSE = 2.72 g C m -2 ;GPP OR2 : R 2 =0.75, RMSE = 2.57 g C m -2 ;GPP OR3 : R 2 =0.82, RMSE =2.20 g C m -2 ). In particular, among these five variants, GPP OR3 The data points produced fit the 1:1 line most closely.

[0096] Step 3-4: Compare the daily changes of the original two-leaf light energy utilization efficiency model and the five variants with the GPP measured by the flux tower during the rice growing season in different years ( Figure 3 ); As can be seen from the figure, the estimated GPP at the ecological site in Jurong, Jiangsu Province is different from the measured GPP EC The results showed similar seasonal changes, with the lowest values ​​occurring at the start and end dates and the peak occurring at the heading stage (DOY: 200-220). OR3 On most days with daily GPP EC The closest, but overestimation is shown late in the 2022 rice growing season. However, GPP AND and GPP OR2 Other models such as showed underestimation in the middle of the rice growing season.

[0097] Step 3-5: Apply the coefficient of determination R 2 , root mean square error RMSE, and relative root mean square error RRMSE to evaluate and verify the results.

[0098] Step 4: Use multi-model process sensitivity analysis to analyze GPP AND 、GPP OR1 、GPP OR2 、GPP OR3 and GPP OR4 Sensitivity analysis was performed on five models and model parameters to assess the impact of model inputs on the uncertainty of GPP estimates and to provide sensitivity measures to process options and model processes. The specific steps included:

[0099] Step 4-1: Use multi-model process sensitivity analysis to analyze GPPAND 、GPP OR1 、GPP OR2 、GPP OR3 and GPP OR4 Sensitivity analysis was performed on five models ( Figure 4 );

[0100] Step 4-2: Use multi-model process sensitivity analysis method to analyze model parameters (APAR sh ,APAR su ,g(CO2),g(T a ) and g(VPD)) to conduct global sensitivity analysis ( Figure 4 ).

[0101] Step 4-3: Sobol' sensitivity analysis of the model input parameters showed that the most sensitive variable among the five input parameters was T a , followed by APAP sh , which can limit rice photosynthesis and regulate the contribution of diffuse radiation to GPP. In addition, the most sensitive process option is GPP OR3 The sensitivity analysis of the two processes showed that the nonlinear effects of environmental factors made a greater contribution to the GPP estimation.

[0102] Step 5: Based on the optimal rice GPP estimation model determined in Steps 3 and 4, for the plot to be estimated, obtain the optimal rice GPP estimation model with the best multi-temporal Sentinel-2 satellite image input to obtain the high-resolution rice GPP and the cumulative GPP during the growing season; specifically including:

[0103] Step 5-1: Determine the optimal model for GPP estimation and the most sensitive model structure and variables based on steps 3 and 4;

[0104] Step 5-2: Determine the optimal GPP estimation model as the improved two-leaf light energy utilization efficiency model mTL-LUE;

[0105] Step 5-3: Use the mTL-LUE model to estimate the cumulative total primary productivity (AGPP) of rice growing season in Jiangsu Province in multiple years, and draw the spatial distribution map of Jiangsu Province for each year (e.g. Figure 5 AGPP during the rice growing season was mainly concentrated in the range of 1600-2400 g C m in 2019. -2 , 1200-2000g C m in both 2020 and 2021 -2 , 1600-2000 g C m in 2022 -2Obvious interannual changes can be seen from the AGPP of rice from 2019 to 2022. The AGPP in 2019 and 2022 is significantly higher than that in 2020 and 2021. In Jiangsu Province, areas with high AGPP values ​​are mainly located in central Jiangsu, such as Xinghua, while areas with low AGPP values ​​are located in southern Jiangsu.

[0106] In a more preferred embodiment, the method further comprises:

[0107] Step 6: Field-level rice yield estimation and accuracy evaluation in Jiangsu Province from 2019 to 2022. The specific steps include:

[0108] Step 6-1: The high spatial resolution GPP accumulation during the rice growing season obtained in step 4 is converted into rice aboveground biomass through carbon utilization efficiency, carbon conversion to biomass coefficient and aboveground biomass allocation coefficient;

[0109] Step 6-2: The converted aboveground biomass is coupled with the PDR-TCT method to estimate the spatiotemporal dynamic index, forming a crop yield estimation method HIDYM driven by high spatial resolution GPP and dynamic harvest index;

[0110] Step 6-3: Conduct a systematic evaluation of HIDYM performance using multiple years of field-measured yield data ( Figure 6 ), it can be seen from the figure that the rice yield data of 4 years are summarized and the dynamic harvest index strategy (R 2 =0.73) is also more accurate than the fixed harvest index (R 2 =0.48), and the accuracy is higher than the single-year yield estimation accuracy. As can be seen from the figure, compared with the fixed harvest strategy, the dynamic harvest index strategy shows stronger convergence on the 1:1 fitting line.

[0111] Step 6-4: Estimation of GPP products with widely used satellite remote sensing GPP MOD17 Compare strategies ( Figure 7 ); Based on GPP OR3 The model's yield estimates were better than those based on GPP in all years. MOD17 Model. Based on GPP MOD17 The model's yield estimate R 2 No more than 0.50 and lower than the GPP-based OR3 The estimated accuracy of the model (R 2 =0.64~0.72). From the GPP MOD17 Model to GPP-based OR3 Model, 2021 R 2 The largest increase (ΔR 2 =0.41), and R in 2022 2 Minimum increase (ΔR 2=0.24). The year with the largest decrease in RMSE for yield estimates was 2021 (▽RMSE=0.47t / ha), while the year with the smallest decrease was 2020 (▽RMSE=0.09t / ha).

[0112] Step 6-5: Use the HIDYM method to estimate rice yield at the field level in Jiangsu Province for multiple years, and generate a multi-year, large-scale, 10-meter-resolution rice yield spatial distribution product such as Figure 8 For the HIDYM method, the estimated rice yield values ​​for the province in 2019 (6.0-15.0 t / ha) and 2022 (6.0-13.0 t / ha) were generally higher than those in 2020 (5.0-12.0 t / ha) and 2021 (5.0-11.0 t / ha). In particular, among the four years, the rice yield in 2019 was the highest, and the rice yield in 2021 was the lowest.

[0113] Step 6-6: Using the coefficient of determination R 2 , root mean square error RMSE, and relative root mean square error RRMSE to evaluate and verify the results.

[0114] The above descriptions are only some embodiments of the present invention. It should be pointed out that a person skilled in the art can make several improvements without departing from the principles of the present invention. These improvements should be regarded as within the protection scope of the present invention.

Claims

1. A method for estimating rice yield at field level by remote sensing based on high-resolution gross primary productivity and dynamic harvest index, comprising the following steps: Step 1: Obtain Sentinel-2 satellite images of the rice growing season in the study area, pre-process and cloud mask the images on the remote sensing cloud platform, and obtain the leaf area index LAI and the start and end dates of the rice growing season; Step 2: Obtain the daily average temperature T during the rice growing season a , relative humidity RH and dew point temperature, saturated water vapor pressure difference VPD; photosynthetically active radiation PAR, atmospheric CO2 data; Step 3: Considering the atmospheric CO2 fertilization effect and environmental factors, determine the optimal model for estimating rice GPP; Step 4: To clarify the environmental factors VPD, T a The physical mechanism of the nonlinear effect of CO2 on GPP AND 、GPP OR1 、GPP OR2 、GPP OR3 and GPP OR4 Conduct sensitivity analysis of the model and model parameters to assess the impact of model inputs on the uncertainty of GPP estimates and provide sensitivity measures to process options and model processes; Step 5: Based on the optimal rice GPP estimation model determined in steps 3 and 4, for the plot to be estimated, obtain the optimal rice GPP estimation model with the best multi-temporal Sentinel-2 satellite image input to obtain the high-resolution rice GPP and the cumulative GPP in the growing season; Step 6: Combine the high-resolution GPP accumulation and dynamic harvest index HI of the rice growing season in step 5 d Coupled with the previous data, the field-level rice yield estimation results were obtained.

2. The method for remote sensing estimation of rice yield at field level based on high-resolution gross primary productivity and dynamic harvest index according to claim 1, characterized in that: Step 1 specifically includes: Step 1-1: Obtain Sentinel-2 images of the rice growing season in the study area on the remote sensing cloud computing platform and perform cloud mask processing; Step 1-2: Use the PROSAIL model and the neural network hybrid model to invert the 8-band reflectance information of the Sentinel-2 image to obtain the rice leaf area index LAI; Step 1-3: Use the rice growing season time series LAI and the inflection point method to calculate the start and end dates of the rice growing season.

3. The method for remote sensing estimation of rice yield at field level based on high-resolution gross primary productivity and dynamic harvest index according to claim 1, characterized in that: Step 2 specifically includes: Step 2-1: Calculate the daily average temperature Ta, relative humidity RH, dew point temperature, and downward shortwave radiation during the entire rice growth period using the start and end dates of the rice growing season calculated in step 1; Step 2-2: Calculate the saturated vapor pressure difference VPD using the daily average temperature, relative humidity and dew point temperature; Step 2-3: Calculate daily photosynthetically active radiation PAR using daily downward shortwave radiation; Steps 2-4: Synchronously acquire atmospheric CO2 data at the National Oceanic and Atmospheric Administration Earth System Research Laboratory.

4. The method for remote sensing estimation of rice yield at field level based on high-resolution gross primary productivity and dynamic harvest index according to claim 1, characterized in that: Step 3 specifically includes: Step 3-1: Use the remote sensing and meteorological data obtained in steps 1 and 2 to drive the original two-leaf light energy utilization model; Step 3-2: Use the bootstrap sampling method to consider the atmospheric CO2 fertilization effect and environmental factors. Use the logic to convert the original two-leaf light energy utilization rate TL-LUE model GPP org Transformed into 5 forms of GPP AND 、GPP OR1 、GPP OR2 、GPP OR3 and GPP OR4 Model, their formulas are as follows; GPP org =(e msu ×APAR su +e msh ×APAR sh )×f(VPD)×f(T amin ) GPP AND =(ε msu ×APAR su +ε msh ×APAR sh )×g(CO2)×g(VPD)×g(T a ) GPP OR1 =(ε msu ×APAR su +ε msh ×APAR sh )×g(CO2)×Min(g(VPD),g(T a )) GPP OR2 =(ε msu ×APAR su +ε msh ×APAR sh )×g(VPD)×Min(g(CO2),g(T a )) GPP OR3 =(ε msu ×APAR su +ε msh ×APAR sh )×g(T a )×Min(g(VPD),g(CO2)) GPP OR4 =(ε msu ×APAR su +ε msh ×APAR sh )×Min(g(CO2),g(VPD),g(T a ))) In the formula, ε msu Indicates the maximum light energy utilization rate of the sun leaf, APAR su Indicates that the sun leaves absorb photosynthetically active radiation, ε msh Indicates the maximum light energy utilization rate of shade leaves, APAR sh represents the photosynthetically active radiation absorbed by the shade leaves, f(VPD) represents the saturated vapor pressure limiting factor, and f(T amin ) represents the minimum temperature limiting factor, g(CO2) represents the carbon dioxide limiting factor, g(VPD) represents the water limiting factor, g(T a ) represents the temperature limiting factor; Step 3-3: Compare the accuracy of the original two-leaf light energy utilization model and five variants with the GPP measured by the flux tower in different years, identify the optimal model for rice GPP estimation, and define the optimal model as the improved two-leaf light energy utilization model mTL-LUE; Step 3-4: Apply the coefficient of determination R 2 , root mean square error RMSE, and relative root mean square error RRMSE were used to evaluate and verify the results, and the optimal model for rice GPP estimation was determined.

5. The method for remote sensing estimation of rice yield at field level based on high-resolution gross primary productivity and dynamic harvest index according to claim 1, characterized in that: Step 4 specifically includes: Step 4-1: Use multi-model process sensitivity analysis to analyze GPP AND 、GPP OR1 、GPP OR2 、GPP OR3 and GPP OR4 Five models were used to conduct model structure sensitivity analysis; Step 4-2: Use the multi-model process sensitivity analysis method to conduct a global sensitivity analysis of the model parameters, including the shade leaf absorption of photosynthetically active radiation APAR sh , the sun leaves absorb photosynthetically active radiation APAR su , carbon dioxide limiting factor g(CO2), temperature limiting factor g(T a ) and water limitation factor g(VPD); Step 4-3: Identify the most sensitive model structure of the five models in step 4-1 and the most sensitive variable in step 4-2.

6. The method for remote sensing estimation of rice yield at field level based on high-resolution gross primary productivity and dynamic harvest index according to claim 1, characterized in that: Step 5 specifically includes: Step 5-1: Determine the optimal model for GPP estimation and the most sensitive model structure and variables based on steps 3 and 4; Step 5-2: Determine the optimal GPP estimation model and the most sensitive model structure as the improved two-leaf light energy utilization efficiency model mTL-LUE; Step 5-3: Use the mTL-LUE model to estimate the cumulative amount of GPP in the local rice growing season for multiple years and draw a local spatial distribution map for each year.

7. The method for remote sensing estimation of rice yield at field level based on high-resolution gross primary productivity and dynamic harvest index according to claim 1, characterized in that: Step 6 specifically includes: Step 6-1: The high spatial resolution GPP accumulation during the rice growing season obtained in step 5 is converted into rice aboveground biomass through carbon utilization efficiency, carbon conversion to biomass coefficient and aboveground biomass allocation coefficient, and the formula is as follows; Among them, B represents the aboveground biomass of rice, GPP t is the GPP of rice from the first day of the growing season to the end of the growing season, is the cumulative gross primary productivity AGPP of rice in the whole growing season; the carbon utilization efficiency CUE adopted is a fixed constant of 0.533, which is used to convert GPP into net primary productivity; F represents the coefficient of carbon conversion into biomass, which is 2.63 for rice; Step 6-2: Convert the spatiotemporal dynamic index HI estimated by the aboveground biomass coupling PDR-TCT method d , forming a crop yield estimation method HIDYM driven by high spatial resolution GPP and dynamic harvest index; Among them, Yield represents rice yield, HI d represents the dynamic harvest index, P is the biomass distribution coefficient of the aboveground part of the plant, which is set to 0.91; W is the grain moisture content coefficient, which is set to 14% for rice.

8. The method for remote sensing estimation of rice yield at field level based on high-resolution gross primary productivity and dynamic harvest index according to claim 1 or 7, characterized in that: Dynamic Harvest Index HI d Obtained by the following steps: Step 601: Based on the multi-temporal Sentinel-2 satellite images within the three key time points of the beginning, peak and end of the rice growing season in the study area, the vegetation index and tassel cap transformation of all Sentinel-2 images within the three key time points are calculated; specifically, the following steps are included: Step 601-1: Obtain information on rice planting distribution and rice growth process in the study area; Step 601-2: Preprocessing and cloud masking of multi-temporal Sentinel-2 images to obtain cloud-free images; Step 601-3: Calculate the Sentinel-2 vegetation index and tasseled cap transformation within the range of the start, peak and end dates of the rice growing season using the rice growth process information in the study area; Step 602: respectively obtaining the vegetation index and tassel cap transformation components of the images at the time points corresponding to the start, peak and end dates of the rice growing season; specifically including: Step 602-1: Extract the start, peak and end dates using the Sentinel-2 vegetation index (EVI) and inflection point method for the study area; Step 602-2: Obtain the Sentinel-2 image vegetation index EVI and LSWI and the three components of tasseled cap transformation, namely brightness, humidity and greenness, for the start, peak and end dates of rice production; Step 603: According to the definition of harvest index and the yield formation process, a phenological difference ratio method PDR is proposed based on Sentinel-2 images of the start, peak and end dates of rice, and then the vegetation index and tasseled cap transformation components are coupled to form a new remote sensing index. A dynamic harvest index estimation model is constructed in combination with the measured harvest index. Finally, the dynamic harvest index estimation model constructed by the vegetation index and tasseled cap transformation components is compared in accuracy by using the exhaustive method and feature screening method to determine the optimal dynamic harvest index estimation model; specifically, the following steps are performed: Step 603-1: Obtain the phenological difference ratio PDR by the following formula: Among them, D SOS , D POS and D EOS vegetation index or tasseled cap transform representing the start, peak, and end dates of the rice growing season, respectively; Step 603-2: coupling the phenological difference ratio PDR with the vegetation index VI and the tasseled cap transform TCT to form a remote sensing index with physiological significance, and then combining it with the measured harvest index to construct a dynamic harvest index estimation model; Step 603-3: using exhaustive method and feature screening method to compare the accuracy of different combinations of vegetation index and tasseled cap transformation components, and determine the optimal dynamic harvest index estimation model; Step 603-4: Apply the coefficient of determination R 2 , root mean square error RMSE, relative root mean square error RRMSE evaluation and verification results; Step 604: For the plot to be estimated, obtain multi-temporal Sentinel-2 satellite images and input the optimal harvest index estimation model to obtain a dynamic harvest index.

9. The method for dynamic harvest index estimation based on multi-temporal high-resolution remote sensing images and phenological information according to claim 8, characterized in that: Vegetation index VI includes enhanced vegetation index EVI and land surface water index LSWI, and the calculation formula is as follows; where ρ NIR , ρ RED , ρ BLUE and ρ SWIR Represent the reflectance values ​​of Sentinel-2 band 8, band 4, band 2 and band 11 respectively; band 8 is 785-900nm, band 4 is 650-680nm, band 2 is 458-523nm, and band 11 is 1565-1655nm; The three components of the Tasseled Cap Transform (TCT) are calculated using the coefficients of the 13 bands of the Sentinel-2 image. The specific calculation formula is as follows: Where m = 1, 2, 3, TCT1 represents the brightness component of the tasseled cap transform TCT, TCT2 represents the greenness component of the tasseled cap transform TCT, and TCT3 represents the humidity component of the tasseled cap transform TCT; B i Represents the reflectance value of the i-th band of the Sentinel-2 image; a i Represents the coefficient of the i-th band.

10. The method for dynamic harvest index estimation based on multi-temporal high-resolution remote sensing images and phenological information according to claim 8, characterized in that: In step 3-2, the remote sensing indicators PDR-TCT and PDR-VI with physiological significance have the following formulas: Among them, TCT SOS 、TCT POS and TCT EOS TCT represents the start date, peak date and end date of the rice growing season m , m = 1, 2, 3, representing the brightness component, green component and humidity component of the tasseled cap transform TCT respectively; Among them, VI SOS VI POS and VI EOS The vegetation indices VI, which represent the start date, peak date, and end date of the rice growing season, include the Enhanced Vegetation Index EVI and the Land Surface Water Index LSWI.

Images