Crop phenology remote sensing extraction method based on shape model local matching

Abstract

The present application provides a kind of crop phenological period remote sensing extraction method based on shape model local matching, comprising: S1, by reference VI time curve and reference phenological period thereon, reference shape model is calculated using SMF method;S2, smooth target VI time curve;S3, reference shape model is matched with the target VI time curve after smoothing processing, and matching function is obtained;S4, phenology estimation: phenological period is estimated by iterative way;S5, the final result is calculated with the target VI time curve after smoothing processing;S6, if correlation coefficient is less than 0.8, then remove phenological period estimation and then output;If correlation coefficient is greater than or equal to 0.8, then directly output.The phenological period is estimated by iterative way, which solves the technical problem that the existing SMF method's phenology estimation, for later phenological period, usually has greater variance than early phenological period, makes the extracted phenological period spatial distribution and time variation produce deviation.

Description

Technical Field

[0001] This invention relates to the field of agricultural remote sensing technology, and in particular to a method for remote sensing extraction of crop phenological periods based on local matching of shape models. Background Technology

[0002] Crop phenology is a critical period in the crop growth process. Accurate crop phenology provides crucial information for agricultural activities such as fertilizer management, irrigation planning, and crop disease prevention. It is also an important parameter for estimating crop yield. Currently, regional-scale crop phenology can only be extracted from satellite remote sensing data. A representative technique for remote sensing extraction of crop phenology is Shape Model Fitting (SMF), proposed by Japanese agricultural remote sensing expert Sakamoto in 2010. The SMF method uses a shape model (h(x)) to perform curve matching between a standard vegetation index (VI) time series with corresponding ground phenological observations and the VI time series of the target pixel (g(x)) through translation and scaling. This is represented as:

[0003]

[0004] In the above formulas, xscale and tshift represent the translation and scaling values, respectively. SMF assumes that the phenological period on the target pixel has the same "relative position" as the phenological period on the matched shape model. Therefore, all phenological periods on g(x) (denoted as ) can be estimated from the corresponding phenological periods on h(x) (denoted as ) using xscale and tshift:

[0005]

[0006] Limitations of SMF Technique: Although the SMF method has been applied in crop phenological stage estimation, it still has two limitations that require the development of new methods to address. First, SMF uses a shape model to match the target VI time series across the entire growing season via linear scaling. This global matching technique assumes that the same xscale and tshift are used for all phenological stages of the crop, which is clearly not the case in reality, as crop growth is not uniform. Furthermore, SMF assumes that the "relative position" of the phenological stages is the same between each target pixel and the matched shape model. Due to this assumption, the time interval ratio between different phenological stages is the same for all pixels. Therefore, SMF cannot be used to analyze relationships between multiple phenological stages. Finally, the scaling factor in the SMF fitting function is related to the phenological stage, which causes the variance of the estimate for a specific phenological stage to depend on that phenological stage itself, which can be expressed as:

[0007]

[0008] Where D(·) and COV(·) represent functions of variance and covariance, respectively.

[0009] In practice, we have found that SMF phenological estimates typically have a larger variance for later phenological periods than for earlier phenological periods, which will cause bias in the spatial distribution and temporal variation of the extracted phenological periods. Summary of the Invention

[0010] To address the shortcomings of existing technologies, this invention provides a remote sensing extraction method for crop phenological periods based on local shape model matching. This method solves the technical problem in existing SMF methods where later phenological periods typically have greater variance than earlier phenological periods, leading to deviations in the spatial distribution and temporal variation of the extracted phenological periods.

[0011] A method for remote sensing extraction of crop phenological periods based on local shape model matching, as described in an embodiment of the present invention, includes:

[0012] S1. Obtain the reference shape model: Ensure that the initial phenological period on the target VI time curve is the same as the reference phenological period on the reference VI time curve, and then use the SMF method to calculate the reference shape model through the reference VI time curve and the reference phenological period on it.

[0013] S2. Smoothing the target VI time curve: The target VI time curve is smoothed using an SG filter;

[0014] S3. Correct the fitting function: Match the reference shape model with the smoothed target VI time curve to obtain the matching function;

[0015] S4. Phenological estimation: The matching function is adapted to a local window, and the phenological period is estimated through iteration to obtain the phenological estimate, and the result remains unchanged until the final result is obtained.

[0016] S5. Output: Output the results of the phenological estimation.

[0017] The technical principle of this invention is as follows: within an adaptive local window, phenological periods are estimated through iteration to obtain phenological estimates, and the final result remains unchanged, thereby making the spatial distribution and temporal variation curves of the estimated phenological periods closer to the actual ones.

[0018] Compared with the prior art, the present invention has the following beneficial effects: by estimating phenological periods through an iterative method, it solves the technical problem in the existing SMF method that the phenological estimation of later phenological periods usually has a larger variance than that of earlier phenological periods, which causes the extracted phenological period spatial distribution and temporal variation to be biased.

[0019] Furthermore, the matching function in step S3 is achieved through:

[0020]

[0021] in, Represents the target VI time curve. Indicates a reference shape model. Indicates the initial phenological period, xscale i tshift represents the translation value. i This indicates the scaling value.

[0022] Furthermore, the phenological estimation in step S4 is calculated using the following formula:

[0023] in, tshift represents all phenological periods on the target VI time curve. i Indicates the scaling value. Indicates the initial phenological period.

[0024] Furthermore, the xscale value in the formula of the matching function... i tshift i The parameters are maximized in a local segment of the target VI time curve. The correlation coefficient between g(x) and g(x) is used to estimate the correlation coefficient, and is expressed as:

[0025]

[0026] Where w represents the half-window width of a local segment of the target VI time curve.

[0027] Furthermore, the half-window width is determined by the local curve characteristics and data noise level of the initial phenological period, and the half-window width is determined by a lookup table method.

[0028] Furthermore, the method for generating the lookup table in the lookup method includes:

[0029] a. Use an SG filter to smooth all VI time curves and aggregate all VI time curves of the same crop within a spatial window;

[0030] b. Randomly add noise to the sample curve and calculate the RMSE values ​​of the noise-simulated VI time curve and the SG-filtered VI time curve in the local interval of the VI time curve.

[0031] c. Estimate the phenological events on the VI time curve by changing the half-window width from 30 to 180 in 15-degree intervals. For a given phenological period, the optimal half-window width value at each noise level is determined as the half-window width value with the smallest difference between the phenological estimate and the reference phenological value of the sample curve, and a lookup table is obtained.

[0032] Furthermore, the noise level in step b is between 5% and 40%, in 5% increments.

[0033] Furthermore, if the correlation coefficient is <0.8, the phenological period estimation in step S4 is removed; if the correlation coefficient is ≥0.8, step S5 is performed directly. Attached Figure Description

[0034] Figure 1 This is a flowchart illustrating the remote sensing extraction method for crop phenological periods based on local shape model matching, according to an embodiment of the present invention.

[0035] Figure 2 This diagram illustrates the impact of stretching operations on the shape model function of the method and the SMF method in this embodiment of the invention.

[0036] Figure 3 This is a standard deviation plot of the phenological estimates of the method and SMF method of this invention in a noise-free simulation scenario for four phenological periods.

[0037] Figure 4 The graph shows the changes in RMSE values, phenological estimation accuracy R² values, and determined half-window width w with noise level for the greening and maturity stages of the method and SMF method in the embodiments of the present invention under a noise simulation scenario.

[0038] Figure 5 This is a spatial distribution map of winter wheat in the North China Plain and the spatial distribution map of ground crop phenology observation stations.

[0039] Figure 6 This is a time series of winter wheat EVI vegetation index and its corresponding nine phenological stages.

[0040] Figure 7 This diagram illustrates the relationship between the field observation phenology and phenological estimation of the method in this embodiment of the invention and the SMF method.

[0041] Figure 8 The variance of phenological periods estimated by the method and SMF method in this embodiment of the invention, and the variance diagram of phenological periods observed at the stations.

[0042] Figure 9 The image shows the phenological results of winter wheat in the North China winter wheat producing area in 2008, plotted using the method and SMF method of this invention.

[0043] Figure 10This is a histogram showing the difference in phenological estimation between the method and the SMF method in this embodiment of the invention, and a time series plot of the average EVI of the pixels at both ends of the histogram.

[0044] The method of the present invention is represented by SMF-S. Detailed Implementation

[0045] The technical solutions of the present invention will be further described below with reference to the accompanying drawings and embodiments.

[0046] like Figure 1 The method for remote sensing extraction of crop phenological periods based on local shape model matching, as shown, includes:

[0047] S1. Obtain the reference shape model: Ensure that the initial phenological period on the target VI time curve is the same as the reference phenological period on the reference VI time curve. The VI time curve is obtained through VI time series data. Then, the SMF method is used to calculate the reference shape model through the reference VI time curve and the reference phenological period on it. The specific initial phenological period on the target VI time curve is extracted through satellite remote sensing data, and the reference phenological period on the reference VI time curve is obtained through actual ground phenological observations.

[0048] S2. Smooth the target VI time curve: Use the SG filter to smooth the target VI time curve.

[0049] S3. Correcting the Fitting Function: Matching the reference shape model with the smoothed target VI time curve to obtain the matching function. The specific matching function is obtained through:

[0050] The formula is obtained; The formula is obtained;

[0051] in, Represents the target VI time curve. Indicates a reference shape model. Indicates the initial phenological period, xscale i tshift represents the translation value. i This indicates the scaling value.

[0052] There are two main differences between the matching function and the fitting function. First, because phenological estimation is mainly controlled by the time dimension, all matching functions remove the parameters yscale and bias. Reducing the number of parameters makes the nonlinear parameter optimization of shape model matching more robust. The nonlinear parameter optimization uses the L-BFGS-B algorithm from the scipy-1.5.0 package in Python 3.8.3. Second, in order to match the initial phenological period... Best to surround The position is stretched and shifted; however, the change in xscale in the fitted function, i.e., xscale i The ×x term caused stretching at the left endpoint, as shown below. Figure 2 As shown, in order to Stretching is applied to both sides of the position, and this is introduced into the fitting function. The term actually compensates for the time shift caused by stretching, such as Figure 2 As shown.

[0053] S4. Phenological Estimation: The matching function is applied using an adaptive local window, and the phenological period is estimated iteratively to obtain the phenological estimate, which remains unchanged until the final result. Specifically, this is achieved through:

[0054] The formula is obtained;

[0055] in, tshift represents all phenological periods on the target VI time curve. i Indicates the scaling value. Indicates the initial phenological period.

[0056] The xscale in the formula of the matching function i tshift i The parameters are maximized in a local segment of the target VI time curve. The correlation coefficient between g(x) and g(x) is used to estimate the correlation coefficient, and is expressed as:

[0057]

[0058] Where w represents the half-window width of a local segment of the target VI time curve.

[0059] The half-window width is determined by the local curve characteristics and data noise level of the initial phenological period. The half-window width is determined by a lookup table method, that is, the initial phenological period is determined by a lookup table method to determine the adaptive local window.

[0060] S5. Output Judgment: If the correlation coefficient R < 0.8, remove the phenological period estimation in step S4; if the correlation coefficient R ≥ 0.8, output directly. This method is used because the noise is too high or the crop type is not correctly identified, and the shape model and the target VI time curve may not be very similar. The correlation coefficient is used to make a judgment.

[0061] The methods for generating the lookup table in the lookup table method include:

[0062] a. Use an SG filter to smooth all VI time curves and aggregate all VI time curves of the same crop within a spatial window.

[0063] b. Randomly add noise to the sample curve, with the specific noise level between 5% and 40% in 5% intervals, and calculate the RMSE value of the VI time curve of the noise simulation and the VI time curve after SG filtering in the local interval of the VI time curve. The larger the RMSE value, the higher the noise level.

[0064] c. Estimate the phenological events on the VI time curve by changing the half-window width from 30 to 180 in 15-degree intervals. For a given phenological period, the optimal half-window width value at each noise level is determined as the half-window width value with the smallest difference between the phenological estimate and the reference phenological value of the sample curve, and a lookup table is obtained.

[0065] The above method yields a lookup table for each phenological stage, allowing w to be determined at any noise level using the lookup table method.

[0066] The method of this invention applies to simulation data, ground phenological observation data, and regional phenological data. Figure 3 Specific tests were conducted in the following aspects:

[0067] Experiment 1: Simulation Experiment Verification

[0068] In the simulation experiments, we used two piecewise logistic functions: the rising phase from the greening period to the VI peak and the falling phase from the VI peak to the dormant period, to simulate the annual VI time series data. Each phase is simulated by a four-parameter logistic function, the formula of which is:

[0069]

[0070] Here, parameters c and d represent the amplitude and background VI value, respectively. represents the date when its amplitude reaches 50%, and b represents the rate of change at . In the simulation, we generated 10,000 annual VI time series data by randomly determining four parameters within a given range. Parameters c and d are in the ranges of [0.5, 0.7] and [0.1, 0.3], respectively. The sum b for the rising phase is randomly determined within the ranges of [80, 120] and [-0.08, -0.05], and the sum b for the falling phase is randomly determined within the ranges of [240, 280] and [0.05, 0.08]. Based on the formula proposed by Shang et al. in 2018, four key phenological periods are defined on the simulated VI time series data, as follows:

[0071]

[0072] Among them, b rising and b descending The parameter b represents the rising and falling phases.

[0073] We compared the method of this invention with the SMF method under noise-free conditions and with different noise levels. Random negative noise of 0%–30% was added to the simulated annual VI time series data. For a fair comparison, both the method of this invention and the SMF method used the same shape model, where each parameter value in the logic function was determined to be the median of the simulation range.

[0074] like Figure 3 As shown in the noise-free simulation experiment results, the method of this invention has a smaller RMSE error value and a higher correlation coefficient between the estimated phenological value and the true value compared with the SMF method.

[0075] For example, the method of this invention estimates the RMSE value for the greening period to be 0.69 days, while the SMF estimates the RMSE value to be 4.42 days. These results indicate that the method of this invention theoretically has a better ability to capture phenological indicators from VI time series data.

[0076] Meanwhile, we calculated the standard deviation (SD) of the phenological estimates using both methods. The results show that the SD of the phenological estimate obtained by the SMF method increases with the increase of the phenological period, underestimating the variance of early phenological periods, such as the greening period. In contrast, the SD of the phenological estimate obtained by the method of this invention is more consistent with the true SD.

[0077] like Figure 4 As shown, the noise simulation experiment results illustrate the changes in the estimated results for the greening and maturity stages relative to the noise level; as... Figure 4 As shown in sections b and e, although the phenological estimation accuracy of both methods decreases with increasing data noise level, the method of this invention outperforms the SMF method at all noise levels, especially at high noise levels. This effect can be attributed to the adaptive local window; to verify this explanation, we tested the method of this invention using a randomly determined half-window width w.

[0078] The results showed that without using an adaptive local window, the phenological estimation error of the method of this invention did indeed increase; at a noise level of 20%, it was even larger than the estimation error of the maturity period by the SMF method. We investigated the half-window width w of the average window size of the method of this invention at different data noise levels. As expected, the half-window width w increased with increasing noise level, although the increasing trend was less pronounced during the greening period. This is because the EVI value is small during the greening period, so the relatively negative noise added to the time series does not significantly reduce the time series, resulting in a small change in the half-window width w during this phenological period.

[0079] Experiment 2: Verification of Ground-based Crop Phenological Observation Data

[0080] In a real-world phenological data experiment, we used the SMF method and the method of this invention to detect the phenological stages of winter wheat in the North China Plain. The distribution of ground crop phenological observation data stations is as follows: Figure 5 As shown, winter wheat accounts for 60% of China's wheat production in this region. Besides winter wheat, other crops include summer maize, soybeans, and rice, which are rotated with winter wheat. Winter wheat in this region is mainly sown after September and harvested before July of the following year, forming a relatively complex annual VI time series with two growth stages, as shown in the figure. Figure 6 As shown, according to field phenological records from the National Meteorological Information Center of the China Meteorological Administration, the nine phenological stages of winter wheat are displayed in [the following text is incomplete and requires further context]. Figure 6 .

[0081] We used the MOD09A1 reflectance product to calculate the EVI (Enhanced Vegetation Index) time series data from 2008 to 2016. Based on the EVI time series data, the winter wheat phenological stages of each station were estimated using the SMF method and the method of this invention. Since field observations only provided the location of the observation stations, to reduce uncertainty, we excluded observation stations where winter wheat pixels accounted for less than 20% of the total pixels within a 20×20 km² area centered on the observation station location. For the remaining field stations, to reduce the spatial scale mismatch between ground-based phenological observations and satellite phenological estimates, we averaged the winter wheat EVI time series data for each station within the 20×20 km² area. The comparison results of the phenological stage estimation results and the ground-based phenological data of the stations are shown below. Figure 7 As shown.

[0082] like Figure 7 As shown, the phenological estimation error of the method of this invention is smaller than the average RMSE value of the SMF method for the nine phenological periods: 9.75 days vs 13.74 days. In particular, the SMF method shows almost no change in the phenological estimates for these early phenological periods, i.e., the variance is underestimated. We further calculated the standard deviation of the phenological estimates, and the results show that the SMF method underestimates the standard deviation SD of the early phenological periods and overestimates the standard deviation SD of the late phenological periods, while the standard deviation SD shows the opposite trend to the actual variation of the phenological periods. Figure 8 As shown, these results demonstrate that the method of the present invention provides higher accuracy in phenological estimation at the site.

[0083] Experiment 3: Verification of Regional Wheat Phenological Period Mapping

[0084] We compared the application of the method of this invention and the SMF method in winter wheat mapping in North China, such as... Figure 9The diagram shows the 2008 germination EMG, spring greening GUD, and maturity date / kernel maturity MD. Through visual comparison, we found that the EMG estimated by the SMF method only varies between 80 and 90 days, while the EMG generated by the method of this invention shows more spatial distribution variation, such as... Figure 9 The a and b parts.

[0085] However, the SMF method observed greater spatial variation for the late phenological period MD; for example, the SMF method estimated the MD from 290 days in the southwest of the region to 315 days in the north, such as... Figure 9 Part c; these spatial observations of phenological changes are consistent with site validation results from SMF that underestimate or overestimate the variance of early or late phenological periods.

[0086] To further compare the two methods, we analyzed the histograms showing the differences between the SMF method and the method of this invention in phenological estimation. The pixels at both ends of the histogram, i.e., the left 10% and the right 10%, are represented as follows: Figure 10 The phenological estimation results shown are averaged and plotted on the EVI time series. We assume that the same phenological period should be in a similar position on the EVI time series.

[0087] The results show that for pixels at both ends of the histogram, the phenological estimates for the same phenological period by SM are not distributed consistently at the same position on the curve. Figure 10 Parts b and c. In contrast, the phenological estimates of the method of the present invention have a more consistent relative position on the EVI vegetation index time series.

[0088] In summary, the method of this invention outperforms the SMF method in all aspects based on the three experimental results.

[0089] Finally, it should be noted that the above embodiments are only used to illustrate the technical solutions of the present invention and are not intended to limit it. Although the present invention has been described in detail with reference to preferred embodiments, those skilled in the art should understand that modifications or equivalent substitutions can be made to the technical solutions of the present invention without departing from the spirit and scope of the technical solutions of the present invention, and all such modifications or substitutions should be covered within the scope of the claims of the present invention.

Claims

1. A remote sensing method for extracting crop phenological periods based on local shape model matching, characterized in that: include: S1. Obtain the reference shape model: Ensure that the initial phenological period on the target VI time curve is the same as the reference phenological period on the reference VI time curve, and then use the SMF method to calculate the reference shape model through the reference VI time curve and the reference phenological period on it. S2. Smoothing the target VI time curve: The target VI time curve is smoothed using an SG filter; S3. Correct the fitting function: Match the reference shape model with the smoothed target VI time curve to obtain the matching function; S4. Phenological estimation: The matching function is adapted to a local window, and the phenological period is estimated through iteration to obtain the phenological estimate, and the result remains unchanged until the final result is obtained. S5. Output: Output the results of the phenological estimation; The matching function in step S3 is obtained through: in, Represents the target VI time curve. Indicates a reference shape model. Indicates the initial phenological period, xscale i tshift represents the translation value. i Indicates the scaling value; The phenological estimation in step S4 is calculated using the following formula: in, This represents all phenological periods on the target VI time curve. Indicates the scaling value. Indicates the initial phenological period.

2. The method for remote sensing extraction of crop phenological periods based on local shape model matching as described in claim 1, characterized in that: The xscale in the formula of the matching function i tshift i The parameters are maximized in a local segment of the target VI time curve. The correlation coefficient between g(x) and g(x) is used to estimate the correlation coefficient, and is expressed as: Where w represents the half-window width of a local segment of the target VI time curve.

3. The method for remote sensing extraction of crop phenological periods based on local matching of shape models as described in claim 2, characterized in that: The half-window width is determined by the local curve characteristics and data noise level of the initial phenological period, and is determined by a lookup table method.

4. The method for remote sensing extraction of crop phenological periods based on local matching of shape models as described in claim 3, characterized in that: The method for generating the lookup table in the table lookup method includes: a. Use an SG filter to smooth all VI time curves and aggregate all VI time curves of the same crop within a spatial window; b. Randomly add noise to the sample curve and calculate the RMSE values ​​of the noise-simulated VI time curve and the SG-filtered VI time curve in the local interval of the VI time curve. c. Estimate the phenology on the VI time curve by changing the half-window width from 30 to 180 in 15 intervals; for a given phenological period, the optimal half-window width value at each noise level is determined as the half-window width value with the smallest difference between the phenological estimate and the reference phenological value of the sample curve, and a lookup table is obtained.

5. The method for remote sensing extraction of crop phenological periods based on local matching of shape models as described in claim 4, characterized in that: The noise level in step b is between 5% and 40%, in 5% increments.

6. The method for remote sensing extraction of crop phenological periods based on local matching of shape models as described in claim 2, characterized in that: If the correlation coefficient is <0.8, the phenological period estimation in step S4 is removed; if the correlation coefficient is ≥0.8, step S5 is performed directly.

Images