A method for remotely estimating rice yield

Abstract

The application provides a rice yield remote sensing estimation method, which comprises the following steps: obtaining NDPI time series of each pixel of a rice plot by using public multispectral remote sensing data; constructing a harvest index HI by using the obtained NDPI time series of each pixel of the rice plot; and establishing a rice yield estimation model by using the constructed harvest index HI and combining histogram information of the maximum NDPI of each pixel in the plot. The rice yield remote sensing estimation method of the application reduces the complexity of the estimation method by constructing a harvest index HI and then establishing the rice yield estimation model of the application, eliminates deviation, improves the application range of the estimation method of the application, and provides a practical technical means for yield estimation of large-area rice plots.

Description

Technical Field

[0001] This method belongs to the field of agricultural remote sensing technology, and in particular relates to a method for estimating rice yield. Background Technology

[0002] Currently, existing methods for estimating rice yield are mainly divided into sampling statistics based on field surveys, methods based on crop models, methods based on empirical models, and methods based on semi-empirical models. For example, the "Review of Data Assimilation of Remote Sensing and Crop Models" uses algorithmic models and specific parameters related to crop growth to simulate the continuous development and growth of crops and estimate yield (see Jin, X., et al., (2018), A review of data assimilation of remote sensing and crop models. European Journal of Agronomy, 92: 141-152, "Review of Data Assimilation of Remote Sensing and Crop Models", European Journal of Agronomy); Marshall and Sharma use spectral bands, spectral indices, backscattering coefficients, radar indices, etc. to construct empirical regression models to estimate yield (see Marshall, M.., et al., (2022), Field-level crop yield estimation with PRISMA and Sentinel-2. ISPRS Journal of Photogrammetry and Remote Sensing, 187: 191-210; Sharma, PK, et al. al., (2022), Assessing the Potentials of Multi-temporal Sentinel-1 SAR Data for Paddy Yield Forecasting Using Artificial Neural Network, Journal of the Indian Society of Remote Sensing, 50(5):895-907.; Lobell used a light use efficiency model to construct a semi-empirical model for yield estimation (see Lobell, DB (2013), The use of satellite data for crop yield gap analysis, Field Crops Research, 143: 56-64).

[0003] However, each of these methods has its own shortcomings: traditional methods obtain rice yield through field surveys or on-site measurements, which not only requires high time and labor costs, but also results in small sample sizes and is prone to bias in the estimation results; crop model-based methods are limited by complex input parameters, and the uncertainty in obtaining canopy state variable data can lead to large errors in estimating crop yield; empirical model-based methods rely on a large number of measured samples for model construction, while ignoring the physiological mechanisms of crops; semi-empirical model-based methods usually require meteorological data for the entire growth period and high temporal resolution remote sensing data as input, and the availability of data greatly limits the application of this type of method. Summary of the Invention

[0004] The technical problem to be solved by the present invention is to provide a remote sensing estimation method for rice yield in order to reduce or avoid the problems mentioned above.

[0005] To address the aforementioned technical problems, this invention proposes a remote sensing estimation method for rice yield, wherein the method includes the following steps:

[0006] First, the NDPI time series of each pixel in the paddy field was obtained using publicly available multispectral remote sensing data;

[0007] Then, using the NDPI time series of each cell in the obtained rice paddy plots, a harvest index is constructed using the following formula. HI :

[0008]

[0009] In the formula, HI represents the harvest index of each pixel; ∑NDPI post This is the cumulative NDPI value after the maximum value at the heading stage. (NDPI value is used as a reference). heading The three NDPI values ​​are used for calculation, representing the cumulative biomass of rice grains during the approximately 30-day maturation period; NDPI is taken as... heading ∑NDPI is calculated using NDPI values ​​from its early growth stages. pre , representing the cumulative biomass of stems and leaves during the growth and reproductive stages of rice;

[0010] Then, the harvest index obtained from the construction was used. HI Based on the histogram information of the maximum NDPI values ​​of each pixel within the plot, a rice yield estimation model is established using the following formula:

[0011]

[0012] In the formula, Yield i HI represents the rice yield of the i-th plot, in kg / mu (unit: kg / mu). iLet be the harvest index of plot i, and be the mean HI of all pixels in that plot; k is the number of intervals in the histogram of the maximum NDPI, and j represents the sequence number for summing each interval, from interval 1 to interval k; NDPI hist_ij Indicates the NDPI in the i-th plot. heading The percentage of pixels whose values ​​fall within the histogram interval, ranging from 0% to 100%; a j Here, b represents the regression coefficient, and b represents the regression intercept term.

[0013] After establishing the rice yield estimation model described above, multiple sample plots within the estimation area were selected to measure their rice yields. The measured rice yields from these sample plots were then input into the estimation model, and the unknown variable 'a' in the model was obtained through least squares regression. j The values ​​of b and the unknown variables are then used to calculate the yield of all rice plots within the estimated area.

[0014] Preferably, the method further includes the following steps: calculating the NDPI curve of each pixel using Sentinel-2 data, and performing cloud removal processing on the NDPI curve using the envelope-constrained curve SG filtering method.

[0015] Preferably, the method further includes the following steps: using Google imagery to segment the fields, and then using the rice index based on synthetic aperture radar data to extract the rice fields from the segmented images.

[0016] Preferably, in the estimation model, k is 4, and the histogram intervals are divided into [0.6, 0.7), [0.7, 0.75), [0.75, 0.8), and [0.8, 1].

[0017] Preferably, the method further includes the step of amplifying rice yield samples using a sample synthesis method:

[0018] The NDPI mean curve of the target plot is synthesized based on the NDPI mean curve of the known sample plots, and the synthesis weight of each sample plot is obtained.

[0019] Weight = [W1…W n ],argmin||X target -Weight·X sample ||1 X target =[NDPI target,1 …NDPI target,t ]

[0020] In the formula, X sample X is a matrix representing the NDPI time series of n known sample plots. target This represents the matrix composed of the NDPI time series of the target plot, where Weight is the composite weight matrix, and W... n NDPI represents the weight of the nth sample plot; t is the time series number. 1,t The NDPI value of the first sample at time t;

[0021] Yield using the known sample plots sample The synthetic yield Yield_snyth of the target plot is calculated using the synthetic weight matrix Weight. target ;

[0022] Yield_snyth target =Weight Yield sample

[0023] In the formula, Yield_snyth target The synthetic yield of the target plot; Weight is the synthetic weight matrix; Yield sample Given the yield matrix of the sample plots, yield n This represents the yield of the nth known sample plot;

[0024] The Yield_snyth of each target plot obtained from the above formula is used. target The synthetic yield represented by this figure serves as the yield in the aforementioned rice yield estimation model. i Substitute them into the estimation model to fit the unknown variables in the estimation model.

[0025] The rice yield remote sensing estimation method of the present invention constructs a harvest index. HI This leads to the establishment of the rice yield estimation model of the present invention, which reduces the complexity of the estimation method, eliminates bias, and expands the application scope of the estimation method of the present invention, providing a practical technical means for yield estimation of large-area rice plots. Attached Figure Description

[0026] The accompanying drawings are intended only to illustrate and explain the present invention and do not limit the scope of the invention.

[0027] Figure 1 The diagram shown is a schematic representation of rice yield estimation results according to a specific embodiment of the present invention.

[0028] Figure 2 What is displayed is Figure 1The rice yield estimation accuracy evaluation curve of the specific embodiment shown. Detailed Implementation

[0029] To provide a clearer understanding of the technical features, objectives, and effects of this invention, specific embodiments are now described. However, those skilled in the art should understand that the following embodiments are not the only limitations on the technical solutions of this invention. Any equivalent transformations or modifications made under the spirit and essence of the technical solutions of this invention should be considered to fall within the protection scope of this invention.

[0030] As described in the background section, traditional methods for surveying or estimating rice yield, based on crop models or empirical models, all suffer from various limitations and biases. Therefore, this invention proposes an improved remote sensing method for estimating rice yield by constructing a harvest index. HI This leads to the establishment of the rice yield estimation model of the present invention, thereby reducing the complexity of the estimation method, eliminating biases, and expanding the application scope of the present invention.

[0031] Specifically, this invention utilizes the Normalized Difference Phenology Index (NDPI) to construct a new harvest index. HI Using the harvest index HI An estimation model for this invention has been established. For the construction and principles of the NDPI index, those skilled in the art can refer to: "Snowless Vegetation Index for Improving the Spring Greening Date of Deciduous Ecosystems," Environmental Remote Sensing (2017), 196, 1-12, Wang Cong et al. The NDPI is an index closely related to vegetation biomass, developed by Wang Cong et al. The time series information of the mean NDPI within a plot represents the combined effect of rice yield on variety and external environmental factors.

[0032] The NDPI time series for each pixel can be calculated using globally available, freely accessible 10-meter resolution Sentinel-2 multispectral remote sensing data. In other words, the estimation method of this invention first requires obtaining the NDPI time series for each pixel of the paddy field using publicly available multispectral remote sensing data.

[0033] To reduce errors, the obtained NDPI time series curves can be preprocessed: NDPI curves for each pixel are calculated using Sentinel-2 data, and the envelope-constrained curve SG filtering method is used to remove clouds from the NDPI curves. For details on the principles and application of SG filtering, please refer to "A Simple Method for Reconstructing High-Quality NDVI Time Series Based on SG Filter," Environmental Remote Sensing, Chen Jin et al., 91, 332-344.

[0034] In addition, to accurately extract paddy fields, Google imagery was used for field segmentation. Then, the segmented images were analyzed using the SAR-based Paddy Rice Index (SPRI) to extract the paddy fields. For details on the principles and calculation of the SPRI index, please refer to "Extracting Robust Indices of Paddy Fields in Cloudy Areas from SAR Time Series," Environmental Remote Sensing, Xu et al., 285, DOI: 10.1016 / j.rse.2022.113374.

[0035] Ultimately, we can obtain NDPI time series data and the extent of rice planting fields.

[0036] Based on the physiological mechanisms of rice, its yield comes partly from the material conversion before the heading stage and partly from the direct photosynthetic products after heading. The NDPI before and after the heading stage can reflect the material storage status of rice during its growth and reproductive stages, thus providing information more relevant to the final rice yield. The cumulative NDPI during the growth period can approximately express the cumulative biomass of stems and leaves before maturity and the cumulative biomass of grains at maturity.

[0037] Based on the above concept, this invention constructs a harvest index using the NDPI time series of each pixel in the obtained rice paddy plot through the following formula. HI :

[0038]

[0039] In the formula, HI represents the harvest index of each pixel; ∑NDPI post This is the cumulative NDPI value after the maximum value at the heading stage. (NDPI value is used as a reference). heading The three subsequent NDPI values ​​are used for calculation, representing the cumulative biomass of rice grains during the approximately 30-day maturation period. The start time of rice growth is determined by the average start time of continuous NDPI increase within the region. Therefore, the NDPI is taken as... heading ∑NDPI is calculated using NDPI values ​​from its early growth stages. pre , representing the cumulative biomass of stems and leaves during the growth and reproductive stages of rice.

[0040] Given the spatial heterogeneity within a plot, each pixel contributes differently to the final yield of that plot. If based on NDPI... heading Histogram distribution information comprehensively considers NDPI at various levels heading The different contributions of pixels to the final average yield of a plot can make full use of the information provided by the spatial heterogeneity of plots, making the yield estimation more accurate.

[0041] Therefore, the present invention further utilizes the above-mentioned harvest index obtained through construction. HI Based on the histogram information of the maximum NDPI values ​​of each pixel within the plot, a rice yield estimation model is established using the following formula:

[0042]

[0043] In the formula, Yield i HI represents the rice yield of the i-th plot, in kg / mu (unit: kg / mu). i Let be the harvest index of plot i, and be the mean HI of all pixels in that plot. k is the number of intervals in the histogram of the maximum NDPI, and j represents the index of each interval from interval 1 to interval k. hist_ij Indicates the NDPI in the i-th plot. heading The percentage of pixels whose values ​​fall within a histogram interval, ranging from 0% to 100%. In a preferred embodiment of the invention, k is 4, and the histogram interval is divided into [0.6, 0.7), [0.7, 0.75), [0.75, 0.8), and [0.8, 1]. j is the regression coefficient, and b represents the regression intercept term.

[0044] After establishing the above rice yield estimation model, multiple sample plots within the estimated area can be selected to measure their rice yields. The measured rice yields from these sample plots can then be input into the estimation model, and the unknown variable 'a' in the model can be obtained through least squares regression. j The values ​​of b and the unknown variables are then used to calculate the yield of all rice plots within the estimated area.

[0045] The larger the survey area of ​​the sample plots, the more model samples are obtained, the more accurate the values ​​of unknown variables are obtained by fitting the model, and the more accurate the estimated rice yield will be. However, as mentioned in the background section, measuring rice yield in the field requires significant time and manpower. Therefore, under the constraints of limited time and manpower, sometimes only a very small number of measured samples can be obtained, making it almost impossible to obtain a reliable and accurate estimation model with such a small sample size. Of course, if there are enough measured samples, this problem does not need to be considered.

[0046] Therefore, for situations where the actual measured sample size is very small, the estimation method of the present invention may further include a step of amplifying the rice yield sample using a sample synthesis method.

[0047] Specifically, the NDPI mean curve of the target plot (sample plots that need expanded production information) can be synthesized based on the known NDPI mean curve of the sample plots to obtain the synthesis weight of each sample plot.

[0048] Weight = [W1…W n ],argmin||X target -Weight·X sample ||1 X target =[NDPI target,1 …NDPI target,t ]

[0049] In the formula, X sample X is a matrix representing the NDPI time series of n known sample plots. target This represents the matrix composed of the NDPI time series of the target plot, where Weight is the composite weight matrix, and W... n Let be the weight of the nth sample plot. t is the time series number, NDPI. 1,t Let be the NDPI value of the first sample at time t.

[0050] Yield using the known sample plots sample The synthetic yield Yield_snyth of the target plot is calculated using the synthetic weight matrix Weight. target .

[0051] Yield_snyth target =Weight Yield sample

[0052] In the formula, Yield_snyth target Yield represents the synthetic yield of the target plot, and Weight is the synthetic weight matrix. sample Given the yield matrix of the sample plots, yield n This represents the yield of the nth known sample plot.

[0053] The Yield_snyth of each target plot obtained from the above formula is used. target The synthetic yield represented by this figure serves as the yield in the aforementioned rice yield estimation model. i Substitute them into the model to fit and estimate the unknown variables in the model.

[0054] The sample synthesis method described above can combine a small number of measured samples from a small number of plots with the target plots by weighting the similarity of their features. This allows for the synthesis of sample yield data for each target plot from a small number of measured samples, thereby expanding the range of samples with comparable accuracy to the measured samples. This method can solve the problem of inaccurate yield estimation over a large area under conditions of insufficient samples, reduce time and labor costs, and achieve higher estimation accuracy.

[0055] The following is a detailed description of an example of rice yield estimation in a part of Wuchang City, Heilongjiang Province, such as... Figure 1 The diagram shown illustrates the rice yield estimation results according to a specific embodiment of the present invention (Wuchang City, Heilongjiang Province). The estimation model provided by the above-described method of the present invention is used to compare the model-fitted yield with the actual sample yield, as shown below. Figure 2 As shown.

[0056] Figure 2 The scatter plot shows the model-fitted yield and the actual sample (including expanded samples) yields based on the land parcel. The scatter points for both the model-fitted yield and the actual sample yield are generally distributed around the 1:1 line. The model's goodness of fit (R²) is... 2 The yield estimation model established in this invention has a relatively reliable accuracy, with a root-mean-square error (RMSE) of 0.61, a mean-absolute error (MAE) of 39.29 kg / mu, and a mean-absolute error (MAE) of 31.67 kg / mu.

[0057] As shown in the table below, the measured yield data of the four plots, without being used in the yield estimation model, show that the yield estimated by the model is very close to the measured yield. Specifically, the root-mean-square error (RMSE) is 28.61 kg / mu, and the mean-absolute error (MAE) is 27.82 kg / mu.

[0058] The table below shows the production data of plots with known plot sizes and the production estimated by the model.

[0059]

[0060] Those skilled in the art should understand that although the present invention has been described with reference to multiple embodiments, not every embodiment contains only one independent technical solution. This description is provided merely for clarity; those skilled in the art should understand the specification as a whole and consider the technical solutions involved in each embodiment as being able to be combined with each other to form different embodiments to understand the scope of protection of the present invention.

[0061] The above description is merely an illustrative embodiment of the present invention and is not intended to limit the scope of the invention. Any equivalent changes, modifications, and combinations made by those skilled in the art without departing from the concept and principles of the present invention should fall within the scope of protection of the present invention.

Claims

1. A remote sensing method for estimating rice yield, characterized in that, The method includes the following steps: First, the NDPI time series of each pixel in the paddy field was obtained using publicly available multispectral remote sensing data; Then, using the NDPI time series of each cell in the obtained rice paddy plots, a harvest index HI is constructed using the following formula: In the formula, HI represents the harvest index of each pixel; ∑NDPI post This is the cumulative NDPI value after the maximum value at the heading stage. (NDPI value is used as a reference). heading The three NDPI values ​​are used for calculation, representing the cumulative biomass of rice grains during the approximately 30-day maturation period; NDPI is taken as... heading ∑NDPI is calculated using NDPI values ​​from its early growth stages. pre , representing the cumulative biomass of stems and leaves during the growth and reproductive stages of rice; Subsequently, using the harvest index HI obtained above, combined with the histogram information of the maximum NDPI of each pixel within the plot, a rice yield estimation model expressed by the following formula was established: In the formula, Yield i HI represents the rice yield of the i-th plot, in kg / mu (unit: kg / mu). i Let be the harvest index of plot i, and be the mean HI of all pixels in that plot; k is the number of intervals in the histogram of the maximum NDPI, and j represents the sequence number for summing each interval, from interval 1 to interval k; NDPI hist_ij Indicates the NDPI in the i-th plot. heading The percentage of pixels whose values ​​fall within the histogram interval, ranging from 0% to 100%; a j Here, b represents the regression coefficient, and b represents the regression intercept term. After establishing the rice yield estimation model described above, multiple sample plots within the estimation area were selected to measure their rice yields. The measured rice yields from these sample plots were then input into the estimation model, and the unknown variable 'a' in the model was obtained through least squares regression. j The values ​​of b and the obtained values ​​of unknown variables are then brought back into the above estimation model, so that the yield of all rice plots in the estimated area can be calculated using this rice yield estimation model. The method further includes the step of amplifying rice yield samples using a sample synthesis method: The NDPI mean curve of the target plot is synthesized based on the NDPI mean curve of the known sample plots, and the synthesis weight of each sample plot is obtained. Weight=[W1 … W n ],arg min||X target -Weight·X sample ||1 X target =[NDPI target，1 … NDPI target，t ] In the formula, X sample X is a matrix representing the NDPI time series of n known sample plots. target This represents the matrix composed of the NDPI time series of the target plot, where Weight is the composite weight matrix, and W... n NDPI represents the weight of the nth sample plot; t is the time series number. 1，t Let be the NDPI value of the first sample at time t; Yield using the known sample plots sample The synthetic yield yield_snyth of the target plot is calculated using the synthetic weight matrix Weight. target ; Yield_snyth target =Weight·Yield sample In the formula, Yield_snyth target The synthetic yield of the target plot; Weight is the synthetic weight matrix; Yield sample Given the yield matrix of the sample plots, yield n This represents the yield of the nth known sample plot; The Yield_snyth of each target plot obtained from the above formula is used. target The synthetic yield represented by this figure serves as the yield in the aforementioned rice yield estimation model. i Substitute them into the estimation model to fit the unknown variables in the estimation model.

2. The method as described in claim 1, characterized in that, The method further includes the following steps: calculating the NDPI curve of each pixel using Sentinel-2 data, and performing cloud removal processing on the NDPI curve using the envelope-constrained curve SG filtering method.

3. The method as described in claim 1, characterized in that, The method further includes the following steps: using Google imagery to segment the fields, and then using the rice index based on synthetic aperture radar data to extract the rice fields from the segmented images.

4. The method as described in claim 1, characterized in that, In the estimation model, k is 4, and the histogram intervals are divided into [0.6, 0.7), [0.7, 0.75), [0.75, 0.8), and [0.8, 1].

Images