A crop phenology detection method in the middle of the season based on stationary satellite data

Abstract

The application 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, acquiring Himawari-8 AHI data, and obtaining a limited VI time sequence through band calculation; performing cloud identification on the AHI limited VI time sequence through a fast cloud detection method to obtain an AHI cloud-removed VI time sequence; performing data synthesis on the AHI cloud-removed VI time sequence through an n-day data 90th percentile synthesis method to generate an AHI n-day VI time sequence, and then smoothing the VI time sequence by using an SG filter; and performing phenological season detection on the smoothed AHI n-day VI time sequence by using an improved SMF-S method. The method can ensure the accuracy of the season detection, greatly reduce the data quantity of the early-stage processing, and is superior to the season detection accuracy of polar-orbit satellite MODIS data with the same spatial resolution.

Description

Technical Field

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

[0002] Crop phenology identifies key physiological stages in crop growth and development, providing crucial information for various agricultural activities such as irrigation planning, fertilization, pest and disease control, and harvest management. It is also an important parameter in crop growth models and crop yield estimation. At regional or global scales, existing studies widely utilize polar-orbiting satellite data for crop phenology detection, commonly including Moderate Resolution Imaging Spectroradiometer (MODIS), Visible Infrared Imaging Radiometer Suite (VIIRS), and Land Remote Sensing Satellite (Landsat). However, frequent cloudy weather conditions can significantly reduce the accuracy of phenological detection using polar-orbiting satellite data. Geostationary satellites provide next-hour observations, which greatly increases the chances of obtaining cloudless observations. Existing studies have demonstrated that geostationary satellite data can be well applied to vegetation phenology monitoring. For example, one study used data from the Spinning Enhanced Visible Infra-Red Imager (SEVIRI) on the Meteosat Second Generation (MSG) geostationary satellite to extract the Normalized Differential Vegetation Index (NDVI) for phenological detection and compared it with the results of MODIS, confirming the potential of geostationary satellite data in phenological research.

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

[0004] Himawari-8 is a new generation of Japanese geostationary weather satellites with 16 observation bands. Its spatial resolution is 0.5 or 1 km in the visible and near-infrared bands, and 2 km in the infrared band. Combined with its shorter imaging time (approximately 10 minutes per scan for the entire satellite and 2.5 minutes per scan for specific regions), many studies have focused on the potential of Himawari-8 AHI data in vegetation phenology monitoring. For example, some studies have found that AHI's high temporal resolution NDVI can better capture seasonal changes in vegetation, potentially improving spring and autumn phenological monitoring and vegetation classification. However, no research has yet addressed the application of AHI data in phenological seasonal detection.

[0005] Most existing vegetation phenology detection methods are post-season detections, conducted after the vegetation growing season ends and a complete VI (vegetation index) time-series curve is obtained. However, agricultural activities require timely and accurate crop phenology information, so in recent years, more and more technicians have begun to study real-time crop phenology detection, also known as mid-season phenology detection. The key issue can be summarized as detecting crop phenology from incomplete VI time-series curves. Liu et al. (2022) proposed the segmented shape model fitting by the Separate phenological stage (SMF-S) method 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, eliminating the dependence of the result variance on phenological periods. It uses an adaptive local window within each phenological period to match the shape model with the target curve, thus ensuring accurate detection even when there are asynchronous changes between different phenological periods. The SMF-S method only requires matching of a local segment of the VI time series, an advantage that suggests its potential for detection during crop phenological seasons. Summary of the Invention

[0006] The purpose of this invention is to overcome the shortcomings of the existing technology and provide a method for detecting crop phenology during the season based on geostationary satellite data, filling the gap in existing research on detecting phenology during the season using geostationary satellite data.

[0007] The technical solution adopted in this invention is as follows:

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

[0009] S1. Determine the study area and time interval, and obtain Himawari-8 AHI data for the current year up to a certain cutoff date; based on the Himawari-8 AHI data, calculate the AHI finite VI time series for the current year through band calculation; at the same time, obtain the spatial distribution of the studied crop and the true values ​​of ground-observed phenology, and generate a shape model.

[0010] S2. Cloud identification is performed on the AHI finite VI time series using a fast cloud detection method. The VI values ​​contaminated by clouds are replaced with null values ​​to obtain the AHI cloud-free VI time series.

[0011] S3. The AHI declouded VI time series obtained in S2 is synthesized by using the 90th percentile of n-day data to generate an AHI n-day VI time series (i.e., the data time interval is n days); the Savitzky-Golay filter is used to correct the negative bias noise of the AHI n-day VI time series to obtain a smoothed VI time series.

[0012] S4. The improved SMF-S method is used to detect the mid-season phenology of the smoothed VI time series in S3, and the mid-season phenology estimate of crops is obtained.

[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 clouds in the AHI finite VI time series, resulting in the AHI de-clouded VI time series; specifically, this includes the following steps:

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

[0015] Six sets of input and output reflectance values ​​are assigned within the range of 0 to 1. Based on these values, cubic spline interpolation is performed to obtain the input-output curve.

[0016] Use the reflectance value of the visible or near-infrared band of AHI data 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. Calculate the Visible-band Cloud Index (VCI) to obtain a preliminary cloud mask;

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

[0019]

[0020] in, , and These represent the enhanced reflectance of AHI bands 1, 3, and 4, respectively. The magnification factor of 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 considered cloudy, thus obtaining a preliminary cloud mask.

[0021] S23. Further identify the missing cirrus cloud pixels in the preliminary cloud mask by using the brightness temperature difference between band 7 and band 13 of each AHI data in the AHI finite VI time series; regard the pixels with a brightness temperature difference between band 7 and band 13 greater than the preset brightness temperature threshold as cloudy, and obtain the final cloud mask.

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

[0023] Furthermore, in step S3, the value of n is 8 in the method for synthesizing the 90th percentile of the n-day data. The sub-diurnal variation of the AHI VI time series is primarily caused by changes in the solar zenith angle and residual cloud noise. Therefore, the AHI de-clouded VI time series obtained in S2 needs further synthesis to obtain a smoother time series. This invention selects the method of synthesizing the 90th percentile of the n-day data to synthesize the original data into an n-day VI time series curve. Experimental verification shows that the method of synthesizing the 90th percentile of data from different zenith angle ranges over 8 days can significantly reduce the amount of data processed in the early stages while maintaining the accuracy of mid-season detection. Then, a Savitzky-Golay filter is used to correct the negative bias noise of the VI time series, smoothing the VI time series.

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

[0025]

[0026] in, The original shape model curve, The shape model curve after scaling and translation. For date, This is the scaling factor. The translation factor is used to calculate the shape model curve after scaling and translation. and The scaling and translation coefficients are obtained from the maximum correlation coefficient of the local segment:

[0027]

[0028]

[0029] in, This represents the Pearson correlation coefficient. It is a predefined time window;

[0030] The final phenological estimate is obtained using the following formula. :

[0031]

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

[0033] This invention presents a method for mid-season phenological detection of crops based on geostationary satellite data. It utilizes Himawari-8 AHI data to detect phenological stages of crops during the growing season based on a finite VI time series, filling a gap in existing research on mid-season phenological detection using geostationary satellite data. First, the method employs a rapid cloud detection method for AHI pixels to obtain a cloud-free AHI VI time series. Then, it uses a 90th percentile synthesis method to synthesize the raw data into an n-day VI time series curve, which is then smoothed using SG filtering. Simultaneously, the method utilizes an improved SMF-S method suitable for mid-season detection scenarios, combined with a proposed Himawari-8 AHI data preprocessing workflow, to provide a phenological mid-season detection scheme suitable for geostationary satellite data. This invention significantly reduces the amount of data processed in the preprocessing stage while maintaining mid-season detection accuracy, and it outperforms the mid-season detection accuracy of MODIS polar-orbiting satellite data with the same spatial resolution. Attached Figure Description

[0034] Figure 1 The spatial distribution map of winter wheat in the North China Plain and the location of ground observation stations (top figure) and the NDVI time series and main phenological stages of winter wheat (bottom figure) are shown in this embodiment of the invention.

[0035] Figure 2 For the Himawari-8 AHI data preprocessing and embodiments of the present invention Methodology flowchart for mid-season phenological detection;

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

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

[0038] Figure 5 AHI and MODIS time series of this invention Scattered map of phenological observation stations during the season;

[0039] Figure 6 A comparison of curves synthesized by time period, synthesized by zenith angle, and synthesized by the 90th percentile of all data in an embodiment of the present invention;

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

[0041] Figure 8 The AHI time series curves synthesized by time period and by zenith angle in this embodiment of the invention are examples of curves synthesized from AHI time series curves. Methods for detecting phenological phenomena during the season. Detailed Implementation

[0042] To make the objectives and technical solutions of the embodiments of this application clearer, the technical solutions of the embodiments of this application will be clearly and completely described below with reference to the accompanying drawings.

[0043] Taking winter wheat from the North China Plain as an example, in this embodiment, Himawari-8 AHI data preprocessing and The process of phenological seasonal detection is as follows: Figure 2 As shown, it includes the following steps:

[0044] S1. Determine the study area and time interval, and obtain Himawari-8 AHI data for the current year up to a certain deadline; based on the Himawari-8 AHI data, calculate the AHI finite VI time series for the current year through band calculation; at the same time, obtain the spatial distribution of the studied crop and the true values ​​of ground-observed phenology, and generate a shape model.

[0045] (1) Observation of ground phenology 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 map above shows the geographical location of the North China Plain, encompassing provinces such as Hebei, Shanxi, Henan, and Shaanxi. The main crop in this region is winter wheat, accounting for 60% of China's total wheat production. Other crops include summer maize, soybeans, and rice, which are rotated with winter wheat. In this region, winter wheat is mainly sown after September of the current year and harvested before July of the following year, forming a relatively complex annual VI time series, such as... Figure 1 As shown in the figure below, winter wheat has two distinct growth stages.

[0046] Based on ground phenological observation records provided by the National Meteorological Information Center of the China Meteorological Administration, the examples delineated nine phenological stages for winter wheat, including sowing, emergence, tillering, overwintering, greening, jointing, heading, milk stage, and maturity. Agronomic descriptions are 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-based phenological observation data for the nine phenological periods in 2016 were obtained. Due to the fact that some ground-based observation stations did not continuously observe certain phenological periods, the number of observation data for each phenological period varies, totaling 378 observation data (Table 1). Among them, the five training sites used to generate the shape model are located at... Figure 1 The above image has been labeled.

[0047] Phenological period Agronomic description quantity Sowing period Plant seeds 42 Seedling stage The first leaf grows from the outer layer of the seed, and is about 2.0 cm long. 43 Tillering stage The tip of the first lateral bud protrudes from the leaf sheath, and is approximately 0.5 to 1.0 cm in length. 40 Overwintering period To ensure maximum yield potential, plants need to have sufficient tillering before entering dormancy. 37 Green period Wheat begins to grow again in winter, with the central leaves growing to about 1.0 to 2.0 cm. 33 Propagation period As the internodes elongate, the initial nodes of the stem have become visible. 46 Heading stage The tip of the spikelet protrudes from the sheath of the flag leaf. 47 milk maturity period The grains in the middle of the ear (including naked oats and the tops of oats) have developed to a normal size and are yellowish-green. 45 Maturity More than 80% of the grains are yellow, as are the glumes and stems; only the upper first and second nodes retain a slight green tinge. 45

[0048] Table 1

[0049] (2) Acquiring Himawari-8 AHI data. In this embodiment, five Himawari-8 AHI data patches covering the North China Plain from July 2016 to July 2017 were acquired from NASA’s GeoNEX website. The patch numbers are h47v04, h48v03, h48v04, h49v03 and h49v04. The data level is L1G and includes top-of-atmosphere (TOA) reflectance data. The time range is UTC 00:00-10:00 (Beijing time 08:00-18:00), with one image every 10 minutes. The total data volume is more than 100,000 images.

[0050] The AHI has 16 channels. Band 3 (red light band) has a spatial resolution of 500m, bands 1, 2, and 4 (blue, green, and near-infrared bands) have a spatial resolution of 1km, and the remaining 12 bands (infrared bands) have a spatial resolution of 2km. This example uses one of the VIs, the Normalized Difference Vegetation Index (NDVI), as the research object. To match the resolution of the near-infrared band, the spatial resolution of band 3 is downsampled to 1km using bilinear interpolation and the weighted average of the nearest 2×2 neighboring pixels for NDVI calculation. (); By using nearest neighbor interpolation, the spatial resolution of bands 7 and 13 is resampled to 1km for subsequent cloud removal processing. The specific bands used in AHI and their functions are shown in Table 2.

[0051] Band number Wavelength (μm) Band Name Spatial resolution (km) use 1 0.47 Blue light (B) 1 Calculate the visible cloud index 3 0.64 Red light (R) 0.5 Calculate NDVI and calculate the visible cloud index. 4 0.86 Near-infrared (NIR) 1 Calculate NDVI and calculate the visible cloud index. 7 3.89 Mid-infrared (MIR) 2 Identifying cirrus clouds 13 10.41 Mid-infrared (MIR) 2 Identifying cirrus clouds

[0052] Table 2

[0053] (3) Acquiring MODIS data. To compare the phenological estimation effect of AHI data in mid-season scenarios, this embodiment selected 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 photosynthetic vegetation activity on Earth's land and supports accurate seasonal and interannual monitoring of Earth's land vegetation. MODIS data has been widely used in vegetation phenological detection. By comparing the phenological detection results with MODIS data, the effectiveness of the Himawari-8 AHI data used in this invention in mid-season phenological detection is verified.

[0054] S2. Cloud identification is performed on the AHI finite VI time series using a fast cloud detection method. Null values ​​are used to replace cloud-contaminated VI values, resulting in 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 rapid cloud detection method is used to identify clouds in the AHI finite VI time series, resulting in the AHI cloud-free VI time series. Specifically, this includes the following steps:

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

[0057] Six sets of input and output reflectance values ​​were assigned within the range of 0 to 1. Based on this, cubic spline interpolation was performed to obtain the input-output curves.

[0058] Use the reflectance value of the visible or near-infrared band of AHI data 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 using bands 1, 3, and 4 of each AHI data frame in the AHI finite VI time series, as shown in the following formula:

[0061]

[0062] in, , and These represent the enhanced reflectance of AHI bands 1, 3, and 4, respectively, with a magnification factor of 255 representing the conversion from reflectance (0-1) to brightness (0-255). Based on the reflectance spectrum, cloud pixels have very low VCI values, while land or water pixels have higher VCI values. Therefore, all pixels with VCI values ​​less than a preset VCI threshold are considered cloudy, resulting in a preliminary cloud mask.

[0063] S23. Further identify the missing cirrus cloud pixels in the preliminary cloud mask by using the brightness temperature difference between band 7 and band 13 of each AHI data in the AHI finite VI time series; regard the pixels with a brightness temperature difference between band 7 and band 13 greater than the preset brightness temperature threshold as cloudy, and 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. The AHI declouded VI time series obtained in S2 is synthesized by using the 90th percentile of n-day data to generate an AHI n-day VI time series (i.e., the data time interval is n days); the Savitzky-Golay filter is used to correct the negative bias noise of the AHI n-day VI time series to obtain a smoothed VI time series.

[0066] Specifically, in the method for synthesizing the 90th percentage of n-day data, n is taken as 8. The sub-diurnal variation of the AHI VI time series is primarily caused by changes in the solar zenith angle and residual cloud noise. Therefore, the AHI de-clouded VI time series obtained in S2 needs further synthesis to obtain a smoother time series. This invention selects the method of synthesizing the 90th percentage of n-day data to synthesize the original data into an n-day VI time series curve. Experiments have shown that the method of synthesizing the 90th percentage of data from 8 days of different zenith angle ranges can significantly reduce the amount of data processed in the early stages while maintaining the accuracy of mid-season detection. Then, a Savitzky-Golay filter is used to correct the negative bias noise of the VI time series, smoothing the VI time series.

[0067] Since agricultural meteorological stations do not provide specific geographical locations of winter wheat fields, this embodiment follows the method of Liu et al. (2022) to reduce the uncertainty of spatial inconsistencies between ground and satellite observations. The steps are as follows: First, winter wheat pixels for 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, to compare with winter wheat phenological observations, the average NDVI time series data of winter wheat pixels within a 20 km × 20 km radius around each station was calculated. Finally, only ground stations with a percentage of winter wheat pixels greater than 20% within the 20 km × 20 km spatial area were retained. The effectiveness of this approach has been demonstrated in previous studies through two experiments: In the first experiment, winter wheat phenological estimates were compared using VI time-series curves averaged at different spatial scales (20 km × 20 km and 5 km × 5 km), revealing very small differences (<2 days) between the estimates; in the second experiment, pixel-level standard deviations were calculated for winter wheat phenological estimates within a 20 km × 20 km spatial area around each station, showing that the standard deviations were much smaller than the phenological estimation errors. These two experiments demonstrate that comparing winter wheat phenology obtained from satellite data with ground-based observations at these stations is acceptable.

[0068] S4. The improved SMF-S method is used to detect the mid-season phenology of the smoothed VI time series in S3, and the mid-season phenology estimate of crops is obtained.

[0069] Specifically, this invention modifies the SMF-S method to adapt to the mid-season phenological detection scenario (SMF-Swithin-season, The shape model generated in step S1 includes shape model curves and phenological periods; for a phenological period on the shape model... Determine the smoothed VI time series of the target pixel. The phenological period For the corresponding local segment, the shape model curve is matched to that local segment using scaling and translation; the scaling and translation formula is as follows:

[0070]

[0071] in, The original shape model curve, The shape model curve after scaling and translation. For date, This is the scaling factor. The translation factor is used to calculate the shape model curve after scaling and translation. and The maximum correlation coefficient of a local segment yields the values ​​of these two coefficients, as shown in the following formula:

[0072]

[0073]

[0074] in, This represents the Pearson correlation coefficient. It is a predefined time window.

[0075] Final phenological estimates We obtain it from the following formula:

[0076]

[0077] The following two experiments are set up to evaluate the effectiveness of the crop phenological seasonal 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 results from AHI vs MODIS phenological mid-season monitoring stations

[0080] To verify whether AHI data can be used for mid-phenological season detection and to assess its accuracy, this example compares the mid-phenological season estimation of NDVI time-series data synthesized from the 90th percentile of all 8 days of AHI data with that of MODIS time-series data. For each phenological period, it is assumed that the NDVI time-series data only captures the phenological period of the shape model. At the location, and obtain different data types. The method yields phenological estimates by calculating the mean absolute error (MAE) between ground observation data and different phenological estimates, and then conducting a quantitative assessment.

[0081] This embodiment uses different data The phenological estimation method was quantitatively evaluated by comparing it with ground observation data provided by winter wheat stations in the North China Plain. Figure 6A scatter plot of all samples between ground-based observation data and phenological estimates is shown. The mean MAE for the nine phenological periods was calculated here; the MAE values ​​for the composite data (90th percentile of all 8-day AHI data) were 11.37 days and 10.85 days, respectively. AHI data performed slightly better in the mid-season detection scenario, and all samples were relatively stable with no outliers (i.e., no significant errors), indicating that the 8-day data synthesized from AHI data using the preprocessing described above can be well used in the mid-season phenological detection scenario.

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

[0083] This embodiment selected the 90th percentile synthesis method to synthesize the raw data into an 8-day NDVI time series curve. To investigate how to select data for synthesis to reduce the processing load of AHI data, the synthesis effects of three different synthesis methods were tested: ① synthesis of all 8-day data using the 90th percentile; ② synthesis of data from different time periods within 8 days using the 90th percentile (08:00-09:00, 09:00-10:00...15:00-16:00); ③ synthesis of data from different zenith angle ranges within 8 days using the 90th percentile (10°-20°, 20°-30°...80°-90°). Station curves were synthesized using these three different synthesis methods. The NDVI time series curves from an average of 5 training stations were used as the shape model curves, and the phenological periods from an average of 5 training stations were used as the shape model phenological periods. Furthermore, the results were analyzed using... The method was used for mid-season phenological detection. 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 presentation compares the curves obtained by combining the 90th percentile of data from different time periods over 8 days (08:00-09:00, 09:00-10:00...15:00-16:00), the 90th percentile of data from different zenith angle ranges over 8 days (10°-20°, 20°-30°...80°-90°), and the 90th percentile of data from all 8 days. A single station's composite curve is shown here. It can be seen that, using the curve obtained by combining the 90th percentile of data from all 8 days as a reference, the curve synthesized based on the zenith angle range is closer in shape, and the shape of the composite curve varies little across different zenith angle ranges, making it more stable.

[0085] Then, the curve similarity R and root mean square error (RMSE) were calculated for the methods of synthesizing the 90th percentage of data from different time periods over 8 days, the methods of synthesizing the 90th percentage of data from different zenith angle ranges over 8 days, and the methods of synthesizing the 90th percentage of data from all data over 8 days. Here, the average value of the results for all stations for each synthesis method was taken, and the results are as follows. Figure 7 As shown, it can be seen that the composite curve based on the 90th percentile of all data over 8 days is generally more similar to the composite curve based on the zenith angle range than the composite curve based on the time range, with a smaller RMSE and a larger R. Among them, the composite curves based on the zenith angle ranges of 10°-20°, 20°-30°, and 30°-40° are even more similar.

[0086] Finally, composite curves of the 90th percentile from data at different time intervals over 8 days and composite curves of the 90th percentile from data at different zenith angle ranges over 8 days were used for analysis. Methods for phenological seasonal estimation, Figure 8 The MAE (Mean Accuracy) of the results is shown. The average accuracy of the detection across nine phenological periods is used. The dashed line in the figure represents the detection accuracy of the curve synthesized using the 90th percentage of all 8 days of data. In comparison, the accuracy of the mid-season phenological estimation results using the time-period synthesis method fluctuates significantly, while the accuracy of the mid-season phenological estimation results using the zenith angle range synthesis method is more stable and generally better than the accuracy of the 90th percentage of all 8 days of data synthesis. Based on the above experimental results, this invention selects a small zenith angle range when synthesizing AHI data, which can significantly reduce the amount of data processed in the early stages while ensuring mid-season detection accuracy.

[0087] 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 method for mid-season detection of crop phenology based on geostationary satellite data, characterized in that, The geostationary satellite data is Himawari-8 AHI data, and the method includes the following steps: S1. Determine the study area and time interval, and obtain Himawari-8 AHI data for the current year up to a certain deadline; based on the Himawari-8 AHI data, calculate the AHI finite VI time series for the current year through band calculation; at the same time, obtain the spatial distribution of the studied crop and the true values ​​of ground-observed phenology, and generate a shape model containing shape model curves and phenological periods. S2. Cloud identification is performed on the AHI finite VI time series using a fast cloud detection method. The VI values ​​contaminated by clouds are replaced with null values ​​to obtain the AHI cloud-free VI time series. S3. The AHI cloud-free VI time series obtained in S2 is synthesized by using the 90th percentile of the n-day data to generate the AHI n-day VI time series; the Savitzky-Golay filter is used to correct the negative bias noise of the AHI n-day VI time series to obtain the smoothed VI time series. S4. The improved SMF-S method is used to detect the mid-season phenology of the smoothed VI time series in S3 and obtain the mid-season estimate of crop phenology. Step S2 includes: S21. Enhance the brightness intensity of the visible and near-infrared bands of each AHI data in the AHI finite VI time series; Six sets of input and output reflectance values ​​are assigned within the range of 0 to 1. Based on these values, cubic spline interpolation is performed to obtain the input-output curve. Use the reflectance value of the visible or near-infrared band of AHI data 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. Calculate the visible cloud index to obtain a preliminary cloud mask; The VCI value of each pixel is calculated using bands 1, 3, and 4 of each AHI data frame in the AHI finite VI time series. The calculation formula is as follows: ； in, , and These represent the enhanced reflectivity of AHI bands 1, 3, and 4, respectively, with a magnification factor of 255 representing the conversion from reflectivity to brightness. All pixels with VCI values ​​less than the preset VCI threshold are considered cloudy to obtain a preliminary cloud mask. S23. In the AHI finite VI time series, the pixels whose brightness temperature difference between band 7 and band 13 of each AHI data frame is greater than the preset brightness temperature threshold are considered as cloudy, and the final cloud mask is obtained. S24. Replace the cloud-contaminated VI values ​​with null values ​​to obtain the AHI cloud-free VI time series; Step S4 includes: The target pixel smoothed VI time series is determined based on the shape model generated in step S1. any phenological period For the corresponding local segment, the shape model curve is matched with that local segment by scaling and translation: ； in, The original shape model curve, The shape model curve after scaling and translation. For date, This is the scaling factor. The translation factor is used to calculate the shape model curve after scaling and translation. and The scaling and translation coefficients are obtained from the maximum correlation coefficient of the local segment: ； ； in, This represents the Pearson correlation coefficient. It is a predefined time window; The final phenological estimate is obtained using the following formula. : 。 2. The method for mid-season crop phenological detection based on geostationary satellite data as described in claim 1, characterized in that, In step S3, the value of n is 8 in the method for synthesizing the 90th percentile of the n-day data.

Images