Crop phenological season detection method based on stationary satellite data

Abstract

The invention discloses a crop phenological season detection method based on stationary satellite data, and belongs to the technical field of remote sensing monitoring. The method comprises the following steps: determining a research area range and a time interval, obtaining Himawari-8AHI data, and obtaining a limited VI time sequence through wave band calculation; carrying out cloud identification on the AHI limited VI time sequence through a rapid cloud detection method to obtain an AHI cloud-removed VI time sequence; performing data synthesis on the AHI cloud removal VI time sequence through a 90-th percentage synthesis method of n-day data to generate an AHI n-day VI time sequence, and smoothing the VI time sequence by using an SG filter; and carrying out phenological season detection on the smoothed AHI n-day VI time sequence by adopting an improved SMF-S method. According to the method, the in-season detection precision can be guaranteed, meanwhile, the data size of earlier-stage processing is greatly reduced, and the in-season detection precision of the polar orbit satellite MODIS data is better than that of the polar orbit satellite MODIS data with the same spatial resolution.

Description

Technical Field

[0001] The invention belongs to the technical field of remote sensing monitoring, and in particular relates to a method for detecting crop phenology in season based on geostationary satellite data. Background Art

[0002] Crop phenology specifies the key physiological stages during crop growth and development and provides key information for various agricultural activities, such as irrigation scheduling, fertilization, pest and disease control, and harvest management. Crop phenology is also an important parameter for crop growth models and crop yield estimation. At the regional or global scale, polar-orbiting satellite data are widely used in existing studies for crop phenology detection, including the Moderate Resolution Imaging Spectroradiometer (MODIS), Visible Infrared Imaging Radiometer Suite (VIIRS), and Land Remote Sensing Satellite (Landsat). However, polar-orbiting satellite data may significantly reduce the accuracy of phenology detection in the case of frequent cloudy weather. Geostationary satellites provide sub-hourly observations, which greatly improves the chances of obtaining cloud-free observations. Existing studies have shown that geostationary satellite data can be well used in vegetation phenology monitoring. For example, one study used data from the Spinning Enhanced Visible Infra-Red Imager (SEVIRI) carried on the Meteosat Second Generation (MSG) geostationary satellite to extract the Normalized Difference Vegetation Index (NDVI) for phenology detection and compared the results with those from MODIS, confirming the potential of geostationary satellite data in phenology research.

[0003] At present, a new generation of geostationary satellite sensors have been launched and put into operation, including Advanced Baseline Imager (ABI) on Geostationary Environmental Satellites (GOES)-16 and -17, Advanced Geosynchronous Radiation Imager (AGRI) on Fengyun (FY)-4A, Advanced Meteorological Imager (AMI) on Geostationary Korea Multi-Purpose Satellite (Geo-KOMPSAT)-2A, Flexible Combined Imager (FCI) on Meteosat Third Generation Imaging (MTG-I) satellite, and Advanced Himawari Imager (AHI) on Himawari-8 and -9. These sensors can image the Earth's hemisphere at intervals of 10-15 minutes, with finer spatial resolutions of 500-1000m, and are equipped with spectral bands suitable for deriving vegetation indices, which can serve as an important source of data for studying vegetation dynamics.

[0004] Himawari-8 is a new generation of Japanese geostationary meteorological satellite with 16 observation bands. Its spatial resolution is 0.5 or 1 km in visible and near-infrared bands, and 2 km in infrared bands. It is equipped with shorter imaging capabilities (about 10 minutes / time for full disk data and about 2.5 minutes / time for regional data). Many studies have focused on the potential of Himawari-8 AHI data in the field of vegetation phenology monitoring. For example, some studies have found that AHI's high temporal resolution NDVI can better capture seasonal changes in vegetation, and has the potential to improve spring and autumn phenology monitoring and improve vegetation classification. However, there is no research on the application of AHI data in phenological season detection.

[0005] Most of the existing vegetation phenology detection is post-season detection, which is carried out after obtaining a complete VI (vegetation index) time series curve after the vegetation growing season ends. Agricultural activities require timely and accurate crop phenology information, so in recent years, more and more technicians have begun to study real-time detection of crop phenology, also known as mid-season phenology detection. The key problem can be attributed to detecting crop phenology from incomplete VI time series curves. Liu et al. (2022) proposed a segmented shape model fitting method (Shape Model Fitting by the Separate phenological stage, SMF-S) based on the Shape Model Fitting (SMF) method (Reference 1: Liu, L., Cao, R., Chen, J., Shen, M., Wang, S., Zhou, J., & He, B. (2022). Detecting crop phenology from vegetation index time-series data by improved shape model fitting in each phenological stage. Remote Sensing of Environment, 277, 113060.). This method modifies the fitting function of the SMF method, eliminates the dependence of the result variance on the phenological period, and uses an adaptive local window in each phenological period to match the shape model with the target curve, thereby ensuring accurate detection even if there are asynchronous changes in different phenological periods. The SMF-S method only requires local segments of the VI time series for matching, and this advantage shows that the SMF-S method has potential in mid-season detection of crop phenology. Summary of the invention

[0006] The purpose of the present invention is to overcome the deficiencies of the prior art and provide a method for mid-season detection of crop phenology based on geostationary satellite data, thereby filling the gap in existing research on mid-season detection of phenology using geostationary satellite data.

[0007] The technical solution adopted by the present invention is as follows:

[0008] A method for detecting crop phenology in season based on geostationary satellite data, characterized in that the geostationary satellite data is Himawari-8AHI data, and the method comprises the following steps:

[0009] S1. Determine the scope and time interval of the study area, obtain the Himawari-8 AHI data of the current year before a certain deadline; based on the Himawari-8 AHI data, obtain the AHI finite VI time series of the current year through band calculation; at the same time, obtain the spatial distribution of the study crops and the ground observation phenological true value to generate the shape model;

[0010] S2, using the fast cloud detection method to identify the cloud in the AHI limited VI time series, replacing the cloud-contaminated VI values ​​with null values, and obtaining the AHI cloud-free VI time series;

[0011] S3, synthesizing the AHI cloud-free VI time series obtained in S2 by the 90th percentile synthesis method of n-day data to generate the AHI n-day VI time series (i.e., the data time interval is n days); using the Savitzky-Golay filter to correct the negative skew noise of the AHI n-day VI time series to obtain a smoothed VI time series;

[0012] S4. Use the improved SMF-S method to perform phenological mid-season detection on the smoothed VI time series in S3 to obtain the mid-season estimation of crop phenology.

[0013] Furthermore, in step S2, since the original AHI data does not provide data quality bands, a fast cloud detection method is used to identify the cloud of the AHI limited VI time series to obtain the AHI cloud-free VI time series; specifically, the following steps are included:

[0014] S21, enhance the brightness intensity of the visible light band and near infrared band of each AHI data in the AHI limited VI time series;

[0015] Assign 6 sets of input reflectance values ​​and output reflectance values ​​in the range between 0 and 1, and perform cubic spline interpolation on this basis to obtain input-output curves;

[0016] The reflectance value of the visible light band or near infrared band of the AHI data is used as input; if the input reflectance value is equal to 0 or 1, the output reflectance value remains unchanged; when the input reflectance value is between 0 and 1, the output reflectance value is the corresponding output value on the input-output curve;

[0017] S22, calculating the visible-band cloud index (VCI) to obtain a preliminary cloud mask;

[0018] The VCI value of each pixel is calculated by band 1, band 3 and band 4 of each AHI data in the AHI finite VI time series. The calculation formula is as follows:

[0019]

[0020] Among them, R1', R3' and R4' represent the enhanced reflectance of AHI band 1, band 3 and band 4 respectively, and the amplification factor 255 represents the conversion from reflectance (0-1) to brightness (0-255); according to the reflectance spectrum, the VCI value of cloud pixels will be very small, while the VCI value of land or water pixels will be higher. Therefore, all pixels with VCI values ​​less than the preset VCI threshold are regarded as cloudy, and the preliminary cloud mask is obtained;

[0021] S23, further identify the cirrus cloud pixels missed in the preliminary cloud mask through the brightness temperature difference between band 7 and band 13 of each AHI data in the AHI finite VI time series; regard the pixels whose brightness temperature difference between band 7 and band 13 is greater than the preset brightness temperature threshold as cloudy, and obtain the final cloud mask.

[0022] S24. Replace the cloud-contaminated VI value with a null value to obtain the AHI cloud-free VI time series.

[0023] Furthermore, in step S3, in the 90th percentile synthesis method of n-day data, the value of n is 8. The sub-daily variation of the AHI VI time series is more caused by the change of the solar zenith angle and the residual cloud noise, so the AHI cloud-free VI time series obtained in S2 needs to be further synthesized to obtain a smoother time series. The present invention selects the 90th percentile synthesis method of n-day data to synthesize the original data into an n-day VI time series curve; through experimental verification, the 90th percentile synthesis method of 8-day data with different zenith angle ranges can greatly reduce the amount of data processed in the early stage while ensuring the mid-season detection accuracy; then the Savitzky-Golay filter is used to correct the negative bias noise of the VI time series and smooth the VI time series.

[0024] Furthermore, in step S4, the SMF-S method was modified to adapt to the mid-season phenology detection scenario (SMF-Swithin-season, SMFS ws ). Specifically, the shape model generated in step S1 includes a shape model curve and a phenological period; a local segment corresponding to any phenological period p0 in the target pixel smoothed VI time series h(x) is determined, and the shape model curve is matched with the local segment by scaling and translating:

[0025]

[0026] Among them, g(.) is the original shape model curve, is the shape model curve after scaling and translation, x is the date, xscale is the scaling factor, and tshift is the translation factor; by calculating the shape model curve after scaling and translation The maximum correlation coefficient with h(x) in the local segment gives the values ​​of the scaling and translation coefficients:

[0027]

[0028] x∈[p0-2×t w ,p0]

[0029] Among them, R(.) represents the Pearson correlation coefficient, t w is a predefined time window;

[0030] The final phenological estimate p'0 is obtained by the following formula:

[0031] p'0=p0-tshift

[0032] Compared with the prior art, the present invention has the following beneficial effects:

[0033] The method of the present invention is a method for mid-season detection of crop phenology based on geostationary satellite data, which is used for Himawari-8AHI data to detect phenological periods of crops within the growing season based on limited VI time series, filling the gap in existing research on mid-season detection of phenology using geostationary satellite data. First, the method of the present invention uses the AHI pixel fast cloud detection method to obtain the AHI cloud-free VI time series; then, the method of the present invention selects the 90th percentile synthesis method to synthesize the original data into an n-day VI time series curve, and smoothes the curve through SG filtering; at the same time, the method of the present invention uses the improved SMF-S method suitable for mid-season detection scenarios, combined with the proposed Himawari-8AHI data preprocessing process, to provide a phenological mid-season detection scheme suitable for geostationary satellite data. The present invention can greatly reduce the amount of data processed in the early stage while ensuring the accuracy of mid-season detection, and is better than the mid-season detection accuracy of polar-orbiting satellite MODIS data with the same spatial resolution. BRIEF DESCRIPTION OF THE DRAWINGS

[0034] Figure 1 The spatial distribution map of winter wheat in the North China Plain and the locations of ground observation stations (upper figure) and the NDVI time series and main phenological periods of winter wheat (lower figure) in an embodiment of the present invention;

[0035] Figure 2 Himawari-8AHI data preprocessing and SMFS of the embodiment of the present invention ws Methods Flow chart of phenological mid-season detection;

[0036] Figure 3 is a reflectivity enhancement curve of an embodiment of the present invention;

[0037] Figure 4The cloud removal effect of the fast cloud detection method of the embodiment of the present invention on the AHI time series of winter wheat pixels;

[0038] Figure 5 SMFS of AHI and MODIS time series in the embodiment of the present invention ws Methods Scatter plot of detection sites during the phenological season;

[0039] Figure 6 The curve comparison of the AHI time series synthesized by time period, synthesized by zenith angle and synthesized by the 90th percentile of all data in the embodiment of the present invention;

[0040] Figure 7 Comparison of the curve similarity R and root mean square error RMSE of the AHI time series synthesized by time period, synthesized by zenith angle and synthesized by the 90th percentile of all data in the embodiment of the present invention;

[0041] Figure 8 The AHI time series of the embodiment of the present invention is synthesized by time period and synthesized by zenith angle SMFS ws Methods Phenological mid-season detection accuracy. DETAILED DESCRIPTION

[0042] In order to make the purpose and technical solution of the embodiments of the present application clearer, the technical solution in the embodiments of the present application will be clearly and completely described below in conjunction with the drawings in the embodiments of the present application.

[0043] Taking winter wheat in the North China Plain as an example, in this example, Himawari-8AHI data preprocessing and SMFS ws Methods The process of phenological mid-season detection is as follows Figure 2 As shown, the following steps are included:

[0044] S1. Determine the scope and time interval of the study area, and obtain the Himawari-8 AHI data of the current year before a certain deadline; based on the Himawari-8 AHI data, obtain the AHI finite VI time series of the current year through band calculation; at the same time, obtain the spatial distribution of the study crops and the true value of ground-based phenology to generate a shape model.

[0045] (1) Ground phenological observation of winter wheat in the North China Plain. The North China Plain is an important geographical region in China, located in northern China (29°-42°N, 105°-122°E). Figure 1The above figure shows the geographical location of the North China Plain, covering Hebei, Shanxi, Henan, Shaanxi and other provinces. The main crop in this region is winter wheat, which accounts for 60% of China's total wheat production. Other crops include summer corn, soybeans and rice, which rotate with winter wheat. In this region, winter wheat is mainly planted after September of the current year and harvested before July of the following year, forming a more complex annual VI time series, such as Figure 1 As shown in the figure below, winter wheat has two distinct growth stages.

[0046] Embodiment According to the ground phenological observation records provided by the National Meteorological Information Center of the China Meteorological Administration, 9 winter wheat phenological periods were delineated, including sowing period, emergence period, tillering period, overwintering period, greening period, jointing period, heading period, milky period and maturity period. The agronomic description is shown in Table 1. Figure 1 The figure below shows the reference positions of the nine phenological periods on the VI time series curve. Ground phenological observation data for the nine phenological periods in 2016 were obtained. Since some ground observation stations did not observe certain phenological periods continuously, the number of observation data for each phenological period was different, totaling 378 observation data (Table 1). The five training stations used to generate the shape model were Figure 1 It is marked in the picture above.

[0047]

[0048] Table 1

[0049] (2) Obtain Himawari-8AHI data. In this embodiment, five blocks of Himawari-8AHI data in the North China Plain from July 2016 to July 2017 were obtained from NASA's GeoNEX website. The blocks are labeled h47v04, h48v03, h48v04, h49v03 and h49v04 respectively; the data level is L1G, including Top-of-Atmosphere (TOA) reflectance data; the time range is UTC 00:00-10:00 (Beijing time 08:00-18:00), one image every 10 minutes, and the total data volume is more than 100,000 images.

[0050] AHI has a total of 16 channels, with a spatial resolution of 500m for band 3 (red light band), 1km for bands 1, 2, and 4 (blue light, green light, and near-infrared bands), and 2km for the remaining 12 bands (infrared bands). The embodiment selected one of the VIs, the Normalized Difference Vegetation Index (NDVI), as the research object. In order to match the resolution of the near-infrared band, the spatial resolution of band 3 is downsampled to 1km by bilinear interpolation and outputting the weighted average of the nearest 2×2 neighboring pixels for the calculation of NDVI. Through nearest neighbor interpolation, the spatial resolution of bands 7 and 13 is resampled to 1 km for subsequent cloud removal. The specific bands used by AHI and their functions are shown in Table 2.

[0051]

[0052] Table 2

[0053] (3) Obtaining MODIS data. In order to compare the phenological estimation effect of AHI data in mid-season scenarios, this embodiment selects MOD13A2 data with the same spatial resolution for comparison. MOD13A2 data provides NDVI bands and data quality bands with a spatial resolution of 1 km and a temporal resolution of 16 days. It can be used to monitor the activities of terrestrial photosynthetic vegetation on Earth and support accurate seasonal and interannual monitoring of terrestrial vegetation on Earth. MODIS data has been widely used in vegetation phenology detection. The phenological detection results of MODIS data are compared to verify the effect of Himawari-8AHI data used in the present invention in phenological mid-season detection.

[0054] S2. Cloud identification is performed on the AHI limited VI time series through a fast cloud detection method, and the cloud-contaminated VI values ​​are replaced with null values ​​to obtain the AHI cloud-free VI time series. (Reference 2: Zhuge, XY, Zou, X., & Wang, Y. (2017). A fast cloud detection algorithm applicable to monitoring and nowcasting of daytime cloud systems. IEEE Transactions on Geoscience and Remote Sensing, 55(11), 6111-6119.)

[0055] Since the original AHI data does not provide data quality bands, a fast cloud detection method is used to identify clouds in the AHI limited VI time series to obtain the AHI cloud-free VI time series. Specifically, the following steps are included:

[0056] S21. Enhance the brightness intensity of the visible and near-infrared bands of each AHI data in the AHI finite VI time series.

[0057] Six groups of input reflectance values ​​and output reflectance values ​​were allocated in the range from 0 to 1, on which basis cubic spline interpolation was performed to obtain the input-output curves.

[0058] The reflectance value of the visible light or near-infrared band of the AHI data is used as input; if the input reflectance value is equal to 0 or 1, the output reflectance value remains unchanged; when the input reflectance value is between 0 and 1, the output reflectance value is the corresponding output value on the input-output curve.

[0059] S22. Calculate the visible-band cloud index (VCI) to obtain a preliminary cloud mask;

[0060] The VCI value of each pixel is calculated by bands 1, 3 and 4 of each AHI data in the AHI finite VI time series. The formula is as follows:

[0061]

[0062] Among them, R'1, R'3 and R'4 represent the enhanced reflectivity of AHI band 1, band 3 and band 4 respectively, and the amplification factor 255 represents the conversion from reflectivity (0-1) to brightness (0-255); according to the reflectivity spectrum, the VCI value of cloud pixels will be very small, while the VCI value of land or water pixels will be higher. Therefore, all pixels with VCI values ​​less than the preset VCI threshold are considered cloudy, and a preliminary cloud mask is obtained.

[0063] S23. Further identify the cirrus cloud pixels missed in the preliminary cloud mask through the brightness temperature difference between band 7 and band 13 of each AHI data in the AHI finite VI time series; the pixels whose brightness temperature difference between band 7 and band 13 is greater than the preset brightness temperature threshold are regarded as cloudy to obtain the final cloud mask.

[0064] S24. Replace the cloud-contaminated VI values ​​with null values ​​to obtain the AHI cloud-free VI time series.

[0065] S3, synthesizing the AHI cloud-free VI time series obtained in S2 by the 90th percentile synthesis method of n-day data to generate the AHI n-day VI time series (i.e., the data time interval is n days); using the Savitzky-Golay filter to correct the negative skew noise of the AHI n-day VI time series to obtain a smoothed VI time series;

[0066] Specifically, in the 90th percentile synthesis method of n-day data, the value of n is 8. The sub-daily variation of the AHI VI time series is more caused by the change of the solar zenith angle and the residual cloud noise. Therefore, the AHI cloud-free VI time series obtained in S2 needs to be further synthesized to obtain a smoother time series. The present invention selects the 90th percentile synthesis method of n-day data to synthesize the original data into an n-day VI time series curve; experiments have proved that the 90th percentile synthesis method of 8-day data with different zenith angle ranges can greatly reduce the amount of data processed in the early stage while ensuring the mid-season detection accuracy. Then, the Savitzky-Gol ay filter is used to correct the negative bias noise of the VI time series and smooth the VI time series.

[0067] Since the agricultural meteorological station does not provide the specific geographical location of the winter wheat field, this embodiment follows the method of Liu et al. (2022) to reduce the uncertainty of the spatial inconsistency between ground and satellite observations. The steps are as follows: First, the winter wheat pixels in 2016-2017 were identified on the North China Plain using the method of Qiu et al. (2017) (Reference 1: Qiu, B., et al., "Winter Wheat Mapping Combining Variations before and after Estimated Heading Dates." ISPRS Journal of Photogrammetry and Remote Sensing 123 (2017): 35-46. Web.); Second, in order to compare with the winter wheat phenological observations, the average value of the NDVI time series data of winter wheat pixels within a 20 km × 20 km range around each station was calculated; Finally, only ground stations with a percentage of winter wheat pixels greater than 20% in the 20 km × 20 km spatial area were retained. The effectiveness of this treatment has been confirmed in two previous studies: in the first experiment, the winter wheat phenology estimates using VI time series curves averaged at different spatial scales (20km×20km and 5km×5km) were compared, and the difference in phenology estimates between them was found to be small (<2 days); in the second experiment, the standard deviation of the winter wheat phenology estimates within the 20km×20km spatial range around each station was calculated at the pixel level, and the standard deviation was found to be much smaller than the phenology estimation error. These two experiments show that it is acceptable to compare the winter wheat phenology obtained from satellite data with the winter wheat phenology observed on the ground at these stations.

[0068] S4. Use the improved SMF-S method to perform phenological mid-season detection on the smoothed VI time series in S3 to obtain the mid-season estimation of crop phenology.

[0069] Specifically, the present invention modifies the SMF-S method to adapt to the mid-season phenology detection scenario (SMF-Swithin-season, SMFS ws ). The shape model generated in step S1 includes a shape model curve and a phenological period; for a phenological period p0 on the shape model, determine the local segment corresponding to the phenological period p0 in the target pixel smoothing VI time series h(x), and match the shape model curve with the local segment by scaling and translation; the scaling and translation formula is as follows:

[0070]

[0071] Among them, g(.) is the original shape model curve, is the shape model curve after scaling and translation, x is the date, xscale is the scaling factor, and tshift is the translation factor; by calculating the shape model curve after scaling and translation The maximum correlation coefficient with h(x) in the local segment obtains the values ​​of these two coefficients, and the formula is as follows:

[0072]

[0073] x∈[p0-2×t w ,p0]

[0074] Among them, R(.) represents the Pearson correlation coefficient, t w is a predefined time window.

[0075] The final phenological estimate p'0 is obtained by:

[0076] p'0=p0-tshift

[0077] The following two experiments are provided to evaluate the effectiveness of the crop phenology mid-season detection method based on geostationary satellite data of the present invention.

[0078] The design and implementation results of each experiment are as follows:

[0079] Experiment 1: Comparison of AHI and MODIS mid-season phenological detection station results

[0080] In order to verify whether AHI data can be used for phenological mid-season detection scenarios and its accuracy, this example compares the NDVI time series data synthesized by the 90th percentile of all AHI 8-day data with the phenological mid-season estimation of MODIS time series data. For each phenological period, it is assumed that the NDVI time series data only obtains the phenological period p0 of the shape model, and obtains SMFS of different data types. ws The phenological estimates of the method were used and the mean absolute error (MAE) between the ground observation data and different phenological estimates was calculated for quantitative evaluation.

[0081] In this embodiment, the SMFS of different data ws The phenological estimates of the proposed method were compared with ground observation data provided by winter wheat stations in the North China Plain for quantitative evaluation. Figure 6The scatter plot of all samples between ground observation data and phenological estimates is shown. The mean of the MAE for the nine phenological periods is calculated here, and the MAE values ​​for the MODIS data and the 90th percentile synthetic data of the AHI 8-day full data are 11.37 days and 10.85 days, respectively. The AHI data performs slightly better in the mid-season detection scenario, and all samples are relatively stable, with no outliers (i.e., no significantly large errors), indicating that the AHI data synthesized for 8 days according to the above preprocessing process can be well used in the phenological mid-season detection scenario.

[0082] Experiment 2: Comparison of different synthesis methods of AHI data

[0083] This embodiment selects the 90th percentile synthesis method to combine the original data into an 8-day NDVI time series curve. In order to study how to select data for synthesis to reduce the processing volume of AHI data, the synthesis effects of three different synthesis methods are tested, namely: ① 90th percentile synthesis of all 8-day data; ② 90th percentile synthesis of data in different time periods of 8 days (08:00-09:00, 09:00-10:00...15:00-16:00); ③ 90th percentile synthesis of data in different zenith angle ranges of 8 days (10°-20°, 20°-30°...80°-90°). The station curves are synthesized by three different synthesis methods, and the NDVI time series curves of 5 training stations are averaged as the shape model curve, and the ground observation phenological periods of 5 training stations are averaged as the shape model phenological period, and then the SMFS is used. ws Methods Mid-season phenology detection was carried out. The results of ②③ were compared with ① through three methods: visual observation, curve similarity comparison, and mid-season detection accuracy evaluation.

[0084] first, Figure 6 The curves of the 90th percentile synthesis method for 8-day data in different time periods (08:00-09:00, 09:00-10:00...15:00-16:00), the 90th percentile synthesis method for 8-day data in different zenith angle ranges (10°-20°, 20°-30°...80°-90°) and the 90th percentile synthesis method for all 8-day data are compared. The synthetic curve of a certain station is selected for display here. It can be seen that with the curve of the 90th percentile synthesis method for all 8-day data as a reference, the shape of the curve synthesized by zenith angle range is closer, and the shape of the synthetic curve for different zenith angle ranges does not change much and is more stable.

[0085] Then, the curve similarity R and root mean square error (RMSE) of the 90th percentile synthesis method of 8-day data in different time periods, the 90th percentile synthesis method of 8-day data in different zenith angle ranges, and the 90th percentile synthesis method of 8-day data were calculated. The average of the results of all station curves for each synthesis method is shown in the following figure. Figure 7 As shown in the figure, it can be seen that taking the curve of the 90th percentile synthesis method of all 8 days of data as a reference, the synthetic curves according to the zenith angle range are generally more similar than the synthetic curves according to the time range, with smaller RMSE and larger R, among which the synthetic curves according to the zenith angle ranges of 10°-20°, 20°-30°, and 30°-40° are more similar.

[0086] Finally, the SMFS was performed using the 90th percentile composite curves of the 8-day data at different time periods and the 90th percentile composite curves of the 8-day data at different zenith angle ranges. ws Method for mid-season estimation of phenology, Figure 8 The result precision MAE is shown. Here, the detection precision of the nine phenological periods is averaged. The dotted line in the figure represents the detection precision of the curve of the 90th percentile synthesis method of all the data for 8 days. In comparison, the precision of the mid-season phenological estimation results synthesized by the time period fluctuates greatly, while the precision of the mid-season phenological estimation results synthesized by the zenith angle range is more stable, and generally better than the precision of the 90th percentile synthesis of all the data for 8 days. Based on the above experimental results, the present invention selects a small range of zenith angles when synthesizing AHI data, which can greatly reduce the amount of data processed in the early stage while ensuring the mid-season detection accuracy.

[0087] Finally, it should be noted that the above embodiments are only used to illustrate the technical solution of the present invention rather than to limit it. Although the present invention has been described in detail with reference to the preferred embodiments, those skilled in the art should understand that the technical solution of the present invention can be modified or replaced by equivalents without departing from the purpose and scope of the technical solution of the present invention, which should be included in the scope of the claims of the present invention.

Claims

1. A method for detecting crop phenology in the middle of the season based on geostationary satellite data, characterized in that: The geostationary satellite data is Himawari-8AHI data, and the method comprises the following steps: S1. Determine the scope and time interval of the study area, obtain the Himawari-8AHI data of the current year before a certain deadline; based on the Himawari-8AHI data, obtain the AHI finite VI time series of the current year through band calculation; at the same time, obtain the spatial distribution of the study crops and the true value of ground-observed phenology, and generate a shape model including shape model curves and phenological periods; S2, using the fast cloud detection method to identify the cloud in the AHI limited VI time series, replacing the cloud-contaminated VI values ​​with null values, and obtaining the AHI cloud-free VI time series; S3, synthesize the AHI cloud-free VI time series obtained in S2 by the 90th percentile synthesis method of n-day data to generate the AHIn day VI time series; use the Savitzky-Golay filter to correct the negative bias noise of the AHIn day VI time series to obtain a smoothed VI time series; S4. Use the improved SMF-S method to perform phenological mid-season detection on the smoothed VI time series in S3 to obtain the mid-season estimation of crop phenology.

2. The method for mid-season crop phenology detection based on geostationary satellite data according to claim 1, characterized in that: Step S2 includes: S21, enhance the brightness intensity of the visible light band and near infrared band of each AHI data in the AHI limited VI time series; Assign 6 sets of input reflectance values ​​and output reflectance values ​​in the range between 0 and 1, and perform cubic spline interpolation on this basis to obtain input-output curves; The reflectance value of the visible light band or near infrared band of the AHI data is used as input; if the input reflectance value is equal to 0 or 1, the output reflectance value remains unchanged; when the input reflectance value is between 0 and 1, the output reflectance value is the corresponding output value on the input-output curve; S22, calculating the visible light cloud index and obtaining a preliminary cloud mask; The VCI value of each pixel is calculated by band 1, band 3 and band 4 of each AHI data in the AHI finite VI time series. The calculation formula is as follows: Among them, R1', R3' and R4' represent the enhanced reflectance of AHI band 1, band 3 and band 4 respectively, and the amplification factor 255 represents the conversion from reflectance to brightness; All pixels with VCI values ​​less than the preset VCI threshold are considered cloudy to obtain a preliminary cloud mask; S23, pixels whose brightness temperature difference between band 7 and band 13 of each AHI data in the AHI finite VI time series is greater than a preset brightness temperature threshold are regarded as cloudy, and a final cloud mask is obtained. S24. Replace the cloud-contaminated VI value with a null value to obtain the AHI cloud-free VI time series.

3. The method for mid-season crop phenology detection based on geostationary satellite data according to claim 2, characterized in that: In step S3, in the 90th percentile synthesis method of n-day data, the value of n is 8.

4. The method for mid-season crop phenology detection based on geostationary satellite data according to claim 3, characterized in that: Step S4 includes: According to the shape model generated in step S1, the local segment corresponding to any phenological period p0 in the target pixel smoothed VI time series h(x) is determined, and the shape model curve is matched with the local segment by scaling and translating: Among them, g(.) is the original shape model curve, is the shape model curve after scaling and translation, x is the date, xscale is the scaling factor, and tshift is the translation factor; by calculating the shape model curve after scaling and translation The maximum correlation coefficient with h(x) in the local segment gives the values ​​of the scaling and translation coefficients: x∈[p0-2×t w ,p0] Among them, R(.) represents the Pearson correlation coefficient, t w is a predefined time window; The final phenological estimate p'0 is obtained by the following formula: p′0=p0-tshift.

Images