Dynamic harvest index estimation method based on multi-temporal high-resolution remote sensing image and phenological information

Abstract

The invention provides a dynamic harvest index estimation method based on a multi-temporal high-resolution remote sensing image and phenological information, and the method comprises the following steps: firstly, obtaining the start date, the vigorous date and the end date of a rice growth season based on a multi-temporal Sentinel-2 image of a remote sensing cloud platform through employing an inflection point method; three components and vegetation indexes of Sentinel-2 image tasseled cap transformation at the three key time points are calculated at the same time; then, a phenological difference ratio method PDR is provided by referring to a harvest index definition, and three components of tasseled cap transformation are coupled with a vegetation index to form a remote sensing index with a physiological basis; and finally, determining an optimal dynamic harvest index estimation model by adopting a feature screening and exhaustion method to realize wide-range multi-year rice harvest index estimation. The method can be used for timely and accurately estimating the harvest index of the large-scale field-level rice, and has great application value in the aspects of accurate estimation of rice yield, rice variety breeding, cultivation technology optimization and effect evaluation.

Description

Technical Field

[0001] The invention belongs to the field of agricultural remote sensing monitoring, and is a dynamic harvest index estimation method based on multi-temporal high-resolution remote sensing images and phenological information. Background Art

[0002] The harvest index (HI) refers to the ratio of grain yield to aboveground biomass at the time of crop harvest. Its essence reflects the distribution ratio of crop assimilation products in grains and nutritional organs. The harvest index is an important agronomic parameter index, which plays an important role in crop yield estimation, variety selection and optimization of production management. Most of the current crop yield estimation studies based on the semi-empirical and semi-mechanistic model of gross primary productivity (GPP) use a fixed value of the crop harvest index. However, as an important indicator that is sensitive to crop types, 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 the model assumes that the proportion of crop aboveground biomass to total biomass and the harvest index are constants, and does not consider the impact of environmental stress on the distribution of crop dry matter between different organs. Therefore, in areas with significant surface heterogeneity or large interannual climate fluctuations, inaccurate harvest index settings will lead to significant deviations in crop yield estimates.

[0003] In order to accurately estimate the harvest index of crops, many scholars have estimated the harvest index based on the ratio of biomass-related variables before and after wheat flowering based on field plot experimental data. In addition, previous researchers have constructed a logarithmic relationship between the harvest index and corn dry matter yield (R 2=0.43), but due to scale issues, this method has not been used in the study of pixel-level harvest index estimation. These methods are theoretically very useful for regional harvest index estimation and mapping using remote sensing information. For example, some scholars used the mean ratio of the normalized difference vegetation index (NDVI) before and after wheat flowering as a characteristic indicator to extract the harvest index of winter wheat, but this method first needs to determine the amplitude range of the harvest index value in the region, which has large uncertainty and high difficulty. Some scholars have proposed a method to estimate the harvest index using NDVI, meteorological data and observed flowering dates. However, these methods face problems such as NDVI saturation effect and pixel-level flowering date estimation. Therefore, it is necessary to propose an accurate, interpretable and efficient remote sensing method to describe the spatiotemporal changes of crop HI in a regional scope. In addition, the extent to which these spatiotemporal differences in HI spatial distribution maps can improve the accuracy of field-level rice yield estimation remains to be solved. Summary of the invention

[0004] The technical problem solved by the present invention is to provide a dynamic harvest index estimation method based on multi-phase high-resolution remote sensing images and phenological information. Through the measured harvest index and agronomic prior knowledge within the province, the three key time points of rice start, vigorous and end dates are determined, the Sentinel-2 image vegetation index and tassel cap transformation are calculated, and a novel phenological difference ratio method PDR is constructed. The three components of vegetation index and tassel cap band are coupled with the PDR method to form a new remote sensing index with a certain physiological basis. Based on the exhaustive method and feature screening method, the optimal spatiotemporal dynamic harvest index estimation model is determined, and a large-scale (provincial) multi-year harvest index estimation and a 10-meter resolution rice harvest index spatial distribution product are realized. This method refers to the definition of the harvest index and has a certain physiological basis. It is a simple, efficient and high-precision dynamic harvest index estimation method, which can be used to draw a spatial distribution map of rice harvest index in different years over a large range.

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

[0006] A dynamic harvest index estimation method based on multi-temporal high-resolution remote sensing images and phenological information comprises the following steps:

[0007] Step 1: Acquire Sentinel-2 satellite images at three key times: the beginning, peak, and end of the rice growing season in the study area. Preprocess and cloud mask the images on the remote sensing cloud platform to calculate the Sentinel-2 vegetation index and tasseled cap transform for the three growth periods. The Sentinel-2 vegetation index includes the enhanced vegetation index (EVI) and the land surface water index (LSWI). The tasseled cap transform includes three components: humidity, greenness, and brightness.

[0008] Step 2: Based on prior knowledge of rice planting and the EVI vegetation index of the multi-temporal Sentinel-2, the inflection point method is used to calculate the start, peak and end dates of the rice growing season, and the vegetation index and tassel cap transformation of the images at the three key time points are obtained respectively;

[0009] Step 3: According to the definition of harvest index and the yield formation process, a phenological difference ratio method (PDR) is proposed based on three Sentinel-2 images of the start, peak and end dates of the rice growing season. Then, the enhanced vegetation index (EVI), the land surface water index (LSWI) and the three components of tasseled cap transformation (humidity, greenness and brightness) are coupled to form a remote sensing index with a certain physiological basis. Then, an estimation model is constructed in combination with the measured harvest index. Finally, the exhaustive method and feature screening method are used to compare the accuracy of the harvest index estimation models constructed by the vegetation index and the tasseled cap transformation components to determine the optimal harvest index estimation model.

[0010] Step 4: For the plot to be estimated, obtain multi-temporal Sentinel-2 satellite images and input the optimal harvest index estimation model to obtain the dynamic harvest index.

[0011] Furthermore, the present invention proposes a dynamic harvest index estimation method based on multi-temporal high-resolution remote sensing images and phenological information, in which the image acquisition time in step 2 is the start date, peak date and end date of rice growth and the phenological difference ratio method defined with reference to the harvest index in step 3.

[0012] Furthermore, the present invention proposes a dynamic harvest index estimation method based on multi-temporal high-resolution remote sensing images and phenological information, wherein the specific steps in step 1 include:

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

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

[0015] Steps 1-3: Use the rice growth process information in the study area to extract the Sentinel-2 vegetation index and tasseled cap transform for the start, peak, and end dates of the rice growing season.

[0016] Furthermore, the present invention proposes a dynamic harvest index estimation method based on multi-temporal high-resolution remote sensing images and phenological information, wherein step 2 specifically comprises the following steps:

[0017] Step 2-1: Use the Sentinel-2 vegetation index EVI and the inflection point method to extract the start date, peak date and end date of rice production in the study area;

[0018] Step 2-2: Use the extracted dates of the three key time points of rice to obtain the Sentinel-2 vegetation index and the three components of the Tasseled Cap Transform: humidity, brightness, and greenness.

[0019] Furthermore, the present invention proposes a dynamic harvest index estimation method based on multi-temporal high-resolution remote sensing images and phenological information, wherein the specific steps in step 3 include:

[0020] Step 3-1: Using the multi-temporal Sentinel-2 images obtained in steps 1 and 2, and referring to the definition of harvest index and the yield formation process, a phenological difference ratio method (PDR) with a physiological basis is proposed;

[0021] Step 3-2: PDR of phenological difference ratio method is coupled with vegetation index and tasseled cap transformation to form remote sensing index with physiological significance, and then combined with measured harvest index to construct dynamic harvest index estimation model;

[0022] Step 3-3: Use the exhaustive method and feature screening method to compare the accuracy of different vegetation indices and tasseled cap transformation components to determine the optimal dynamic harvest index estimation model;

[0023] 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.

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

[0025] 1. The dynamic harvest index estimation method based on multi-temporal high-resolution remote sensing images and phenological information of the present invention refers to the definition of harvest index and yield formation process, and has strong mechanism and universality;

[0026] 2. The dynamic harvest index estimation method based on multi-temporal high-resolution remote sensing images and phenological information of the present invention only requires remote sensing images at three time points, making it an efficient method with low uncertainty and simple implementation;

[0027] 3. The dynamic harvest index estimation method based on multi-temporal high-resolution remote sensing images and phenological information of the present invention can be used as an accurate and interpretable method to effectively capture the differences in rice harvest index in different years in a large area. BRIEF DESCRIPTION OF THE DRAWINGS

[0028] Figure 1 It is the technical roadmap of the dynamic harvest index estimation method based on multi-temporal high-resolution remote sensing images and phenological information of the present invention.

[0029] Figure 2 It is a schematic diagram of key time point extraction of the dynamic harvest index estimation method based on multi-temporal high-resolution remote sensing images and phenological information of the present invention.

[0030] Figure 3 It is a schematic diagram of the construction principle of the phenological difference ratio method of the dynamic harvest index estimation method based on multi-temporal high-resolution remote sensing images and phenological information of the present invention.

[0031] Figure 4 The PDR-TCT method and the remote sensing cumulative value method are the accuracy comparisons of the dynamic harvest index estimation method based on multi-temporal high-resolution remote sensing images and phenological information, where (A) and (B) represent the optimal tassel cap transformation variable combination and PDR of rice in 2020, and the cumulative GPP (AGPP) after POS, respectively. post ) and the cumulative GPP before POS (AGPP pre ) ratio.

[0032] Figure 5 This is the spatial distribution map of the dynamic harvest index of rice in Jiangsu in 2020 estimated by (A) the PDR-TCT indicator method and (B) the environmental stress method.

[0033] Figure 6 This is the spatial distribution map of the dynamic harvest index of rice in Jiangsu Province (AD) from 2019 to 2022 estimated based on the PDR-TCT indicator method. DETAILED DESCRIPTION

[0034] 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.

[0035] The present invention is implemented based on Sentinel-2 satellite images. The specific usage characteristics are shown in Table 1.

[0036] Table 1 Sentnel-2 characteristic variables

[0037]

[0038] 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 training and testing of the field-scale rice harvest index estimation model were carried out in 2020, and the traditional remote sensing cumulative value and environmental stress method were compared, which can effectively verify the advancement, universality and robustness of the present invention.

[0039] A dynamic harvest index estimation method based on multi-temporal high-resolution remote sensing images and phenological information, such as Figure 1 As shown in the figure, it mainly includes data preprocessing, collection of field-scale rice harvest index, determination of the start, peak and end dates of the rice growing season, proposal of the phenological difference ratio method, determination of the optimal remote sensing index of PDR-TCT, comparison and evaluation of the accuracy of the remote sensing cumulative value method and the environmental stress method, and mapping of the spatial distribution of rice harvest index in large areas and years. It specifically includes the following steps:

[0040] Step 1: Acquire Sentinel-2 satellite images at different periods during 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, vegetation index, and tasseled cap transformation calculations. 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 performs cloud removal on the original Sentinel-2 surface reflectance image based on this dataset to obtain a cloud-free Sentinel-2 image; then, the present invention calculates the two vegetation indices EVI and LSWI and the tasseled cap transformation brightness (B), humidity (W), and greenness (G) based on the original Sentinel-2 bands. Specifically:

[0041] The calculation formulas of the two vegetation indices, EVI and LSWI, are as follows;

[0042]

[0043] where ρ NIR , ρ RED , ρ BLUE and ρ SWIR They 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.

[0044] 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:

[0045]

[0046] 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. As shown in Table 2:

[0047] Table 2: Comparison table of coefficient values ​​for 13 bands of Sentinel-2 images

[0048]

[0049]

[0050] Step 2: Determine the start, peak, and end dates of the rice growing season based on Sentinel-2 images within the three key growth periods, including:

[0051] Step 2-1: Use the Sentinel-2 vegetation index EVI and the inflection point method in the study area to extract the start date, peak date, and end date of rice. Specifically, the inflection point method (i.e., the second-order derivative is zero) is used to estimate the start and end dates (DOY1 and DOY2) of the rice growing season in Jiangsu Province, and the root mean square error between the estimated start date and the observed date of the rice growing season in 2020 is 4.9 days, thus drawing the spatial distribution map of the start and end dates of the rice growing season in Jiangsu Province in 2020;

[0052] Step 2-2: Use the extracted dates of the three key time points of rice to obtain the Sentinel-2 vegetation index and the three components of the Tasseled Cap Transform, brightness, humidity, and greenness.

[0053] Step 3: According to the definition of harvest index and the yield formation process, a phenological difference ratio method (PDR) is proposed based on three Sentinel-2 images of the start, peak and end dates of the rice growing season. Then, the enhanced vegetation index (EVI), the land surface water index (LSWI) and the tasseled cap transform (TCT) three components of humidity, greenness and brightness are coupled to form a remote sensing index with a certain physiological basis. Then, an estimation model is constructed in combination with the measured harvest index. Finally, the exhaustive method and feature screening method are used to compare the accuracy of the harvest index estimation model constructed by combining the vegetation index and the tasseled cap transform three components to determine the optimal harvest index estimation model, which includes:

[0054] Step 3-1: The remote sensing images of rice start date, peak date and end date obtained in steps 1 and 2 are used to refer to the harvest index definition, and a phenological difference ratio (PDR) method is proposed. The formula is as follows;

[0055]

[0056] Where D SOS , D POS and D EOS The formula requires remote sensing information of three key growth stages of rice (SOS, POS and EOS) as input data. Figure 3 .

[0057] Step 3-2: PDR method is coupled with three bands, namely, enhanced vegetation index (EVI), land surface water content (LSWI) and tasseled cap transform, to form a new remote sensing index with a certain physiological basis. The formula is as follows. Then, a dynamic harvest index estimation model is constructed in combination with the measured harvest index.

[0058]

[0059] Among them, TCT SOS 、TCT POS and TCT EOS The three components TCT representing the SOS, POS and EOS dates respectively m , m=1, 2, 3, respectively representing the brightness component of the tasseled cap transform TCT, representing the greenness component of the tasseled cap transform TCT, and representing the humidity component of the tasseled cap transform TCT.

[0060]

[0061] Among them, VI SOS VI POS and VI EOS The vegetation indices VI, including the Enhanced Vegetation Index EVI and the Land Surface Water Index LSWI, represent the SOS, POS and EOS dates respectively.

[0062] Step 3-3: The exhaustive method and feature screening method are used to compare the accuracy of the harvest index estimation model constructed by combining the vegetation index and the tasseled cap transformation to determine the optimal harvest index estimation model as shown in Table 3. It can be seen from the table that the combination C1(G+W)(R 2 =0.77, RMSE = 0.020) in estimating rice HI d The performance of the LSWI is better than that of the vegetation index, single TCT variables and other combinations. It is worth noting that the variables related to moisture (LSWI: R 2=0.69, RMSE=0.023; W: R 2 =0.71, RMSE = 0.021) in estimating HI d The performance of TCT was second only to the best combination, and the variables of TCT were better than the vegetation index LSWI. However, the greenness-related variables (EVI: R 2 =0.37, RMSE=0.032; G: R 2 =0.40, RMSE = 0.031) to estimate HI d The accuracy is slightly lower, and the vegetation index EVI is slightly lower than the TCT variable;

[0063] Table 3 Accuracy comparison chart

[0064]

[0065] 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.

[0066] Step 4: For the plot to be estimated, obtain multi-temporal Sentinel-2 satellite images and input the optimal harvest index estimation model to obtain the dynamic harvest index.

[0067] Furthermore, the embodiment compares the remote sensing cumulative value method and the environmental stress method with the proposed method based on the PDR-TCT index, and the specific steps include:

[0068] 1): The remote sensing cumulative value used by previous researchers is 1-NDVI post / NDVI pre Index transformation to obtain 1-AGPP post / AGPP pre Combined with the measured harvest index, a dynamic harvest index estimation model is constructed. Figure 4 , the PDR-TCT indicator-based method was used in the 2020 rice HI d (R 2 =0.43, RMSE = 0.031) outperforms 1-AGPP post / AGPP pre ;

[0069] 2): Using the environmental stress method proposed by previous researchers, the spatial distribution map of rice harvest index in Jiangsu Province in 2020 was drawn as follows Figure 5 , we can see that the HI of the former d The spatial distribution map shows that Jiangsu Province has obvious spatial variation in the north-south direction, and the HI d The spatial distribution map shows a relatively uniform distribution across the province, with most values ​​concentrated in the range of 0.40 to 0.45;

[0070] 3): The two methods of step 1) and step 2) were compared with the method based on the PDR-TCT indicator to better reflect the advantages of the PDR-TCT method;

[0071] 4): Application determination coefficient R 2 , root mean square error RMSE, relative root mean square error RRMSE evaluation and verification results;

[0072] 5): The dynamic harvest index model based on the PDR-TCT indicator method was used to realize the spatial distribution map of the harvest index in Jiangsu Province from 2019 to 2022. Figure 6 From the figure, we can see that rice HI d In 2019, the main range was 0.39-0.56, in 2020 and 2021 it was mainly 0.31-0.51, and in 2022 it was 0.40-0.52. The rice HI from 2019 to 2022 showed obvious interannual changes, and the HI in 2019 and 2022 was significantly higher than that in 2020 and 2021. In Jiangsu Province, there was no obvious pattern in the high and low HI areas.

[0073] 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 dynamic harvest index estimation method based on multi-temporal high-resolution remote sensing images and phenological information, comprising the following steps: Step 1: 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 transform of all Sentinel-2 images within the three key time points are calculated; Step 2: Obtain 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 respectively; Step 3: 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. Then, the vegetation index and tasseled cap transformation 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. Step 4: For the plot to be estimated, obtain multi-temporal Sentinel-2 satellite images and input the optimal harvest index estimation model to obtain the dynamic harvest index.

2. The method for dynamic harvest index estimation based on multi-temporal high-resolution remote sensing images and phenological information according to claim 1, characterized in that: Step 1 specifically includes: Step 1-1: Obtain information on rice planting distribution and rice growth process in the study area; Step 1-2: Preprocessing and cloud masking of multi-temporal Sentinel-2 images to obtain cloud-free images; Steps 1-3: Use the rice growth process information in the study area to calculate the Sentinel-2 vegetation index and tasseled cap transform within the start, peak, and end dates of the rice growing season.

3. The method for dynamic harvest index estimation based on multi-temporal high-resolution remote sensing images and phenological information according to claim 2, 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 They 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.

4. The method for dynamic harvest index estimation based on multi-temporal high-resolution remote sensing images and phenological information according to claim 2, characterized in that: 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.

5. The method for dynamic harvest index estimation based on multi-temporal high-resolution remote sensing images and phenological information according to claim 1, characterized in that: Step 2 specifically includes: Step 2-1: Extract the start, peak and end dates using the Sentinel-2 vegetation index (EVI) and inflection point method for the study area; Step 2-2: Obtain the Sentinel-2 image vegetation index EVI and LSWI as well as the three components of tasseled cap transformation (brightness, humidity and greenness) for the start, peak and end dates of rice production.

6. The method for dynamic harvest index estimation based on multi-temporal high-resolution remote sensing images and phenological information according to claim 1, characterized in that: Step 3 specifically includes: Step 3-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 3-2: Couple 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 combine it with the measured harvest index to construct a dynamic harvest index estimation model; Step 3-3: Use the exhaustive method and feature screening method to compare the accuracy of the combination of different vegetation indices and the three components of the tasseled cap transformation to determine the optimal dynamic harvest index estimation model; 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.

7. The method for dynamic harvest index estimation based on multi-temporal high-resolution remote sensing images and phenological information according to claim 1, 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