General high-resolution rice field planting intensity and calendar surveying and mapping system

Abstract

A system for universal and automatic high resolution rice planting intensity and calendar mapping employs one or more processors for performing the steps of acquiring remote C-band Synthetic Aperture Radar (SAR) data to extract potential transplanting signals of a rice field. Optical, infrared and passive microwave remote sensing data are also collected and used for checking the authenticity of potential transplanting signals. A map of rice planting intensity and planting calendar (rice transplanting date and harvesting date) is generated based on authenticated transplanting signals without using ground live data, such as field reference samples and statistical data.

Description

Technical Field

[0001] The present invention relates to a system for mapping rice paddy crop intensity and the planting and harvesting dates of paddy rice, and more specifically, to automatically performing the mapping at a relatively high resolution. Background Technology

[0002] The United Nations aims to eliminate hunger, ensure food security, and promote sustainable agriculture worldwide by 2030. Rice is an important staple food, especially for Asians. Given the significant uncertainty in statistical data for most developing countries (Amalia and Kadir, 2019), consistent spatial and temporal monitoring of rice harvested area is crucial for future food security and agricultural output in major rice-producing countries, particularly those with reduced labor or those sacrificing paddy fields to reduce net carbon emissions (Wang et al., 2017; Zhang et al., 2017). The “rice planting calendar” primarily refers to the number of days in a year for planting and harvesting rice (DOY) and can significantly influence rice yield (Ding et al., 2020; Hashemi et al., 2022). Furthermore, high-quality data on rice harvest area and planting calendar are crucial for improving estimates of irrigation water demand and farmland greenhouse gas emissions (e.g., methane (CH4) emissions), which are key components of water conservation and climate change mitigation (Carlson et al., 2017; Wang et al., 2022; Zhang et al., 2020).

[0003] Since most paddy rice fields are flooded before crop emergence, they can be separated from dry land (i.e., non-paddy fields) by referencing multi-temporal remote sensing imagery. The flooding period (i.e., transplanting period) of rice is typically identified by combining optically based moisture and vegetation indices, with Moderate Resolution Imaging Spectroradiometer (MODIS) data being the most commonly used (Dong and Xiao, 2016; Han et al., 2022; Murali Krishna et al., 2011; Zhang et al., 2017). However, most farmland in Massachusetts has a field size of less than 6400 m². 2(80m × 80m), and some paddy fields are even less than 10m wide (Shen et al., 2023), far below the spatial resolution of the relevant MODIS bands (500 to 1000m). The mixed pixel problem can lead to significant uncertainty in MODIS-based rice mapping results (Liu et al., 2019). Although higher-resolution optical images (e.g., Landsat, Sentinel-2) are also available, they are more susceptible to interference from clouds, shadows, and other noise due to their lower temporal frequency (5–16 days), especially in areas where rice flooding coincides with the rainy season (Singha et al., 2019). Therefore, high-resolution rice mapping using Landsat and / or Sentinel-2 is only robust in tropical / subtropical monsoon Asia (Wei et al., 2022).

[0004] Due to the strong penetrating power of microwaves and their near-independence from weather conditions, an increasing number of studies are incorporating microwave data into rice mapping in cloudy areas (Shen et al., 2023; Song et al., 2018). Active microwave data, especially synthetic aperture radar (SAR) data, plays a crucial role in high-resolution rice mapping. Fully submerged crops have smooth surfaces, which facilitate specular reflection of radar signals (Guan et al., 2023). Therefore, the SAR backscattering coefficient (backscatterer) of paddy fields during transplanting is typically significantly lower than that of other land surfaces (e.g., soil, vegetation, impermeable land) and close to or slightly higher than that of water bodies (Xu et al., 2023). Thresholding methods (Han et al., 2021), exponential construction methods (Xu et al., 2023), time-series similarity-based methods (Pan et al., 2021), and machine learning methods (Sun et al., 2023) have been applied to SAR data for rice mapping. These methods ignore the spatial heterogeneity of SAR backscattering thresholds and the impact of soil moisture changes or freezing on SAR backscattering, heavily relying on regional statistics and ground-based data. These drawbacks limit their applicability to specific regions rather than continental or global scales (Shen et al., 2023; Tian et al., 2023). The lack of a universal automatic mapping algorithm makes it impossible to obtain high-resolution (<30m) rice maps with consistent spatial and temporal coverage globally.

[0005] The well-known global rice planting calendar, RiceAtlas, is primarily compiled using regional statistical data and is therefore provided at the national / provincial / regional level (Laborte et al., 2017). Recent studies have utilized remote sensing data, including MODIS, Sentinel-2, and SAR data, to derive planting and harvesting dates for rice within each grid. However, due to the lack of reliable high-resolution rice distribution maps and the relatively high uncertainty in their grid-scale rice calendar estimates, these studies have had to aggregate rice calendar data to the provincial / regional level or with very coarse resolution (e.g., 0.5°) (Li et al., 2023; Mishra et al., 2021; Zhang et al., 2022; Zhao et al., 2023b). Considering the high spatial heterogeneity of rice calendars in countries such as India (More et al., 2016), a high-resolution gridded rice calendar map covering the entire world is essential.

[0006] Chinese patent document CN110472184B discloses a system for mapping rice using only optical remote sensing data (Landsat multispectral imagery). Therefore, it can only be applied to specific areas where cloud cover is not very frequent. This method uses field reference samples to calculate a local threshold for the rice photometric index, which is unsuitable for areas without sufficient field reference values. Chinese patent document CN109345555B discloses a system for mapping rice using commercial SAR data (COSMO-SkyMed imagery). It acquires field reference samples from commercial high-resolution satellite imagery and uses these reference samples as training data for supervised classification. It would be beneficial to achieve accurate mapping without the need to acquire field reference samples. Summary of the Invention

[0007] This invention relates to a general and automated high-resolution rice planting intensity and calendar mapping system that uses C-band SAR data to extract potential rice paddy planting signals and verifies the authenticity of these signals using optical, infrared, and passive microwave remote sensing data. In particular, it relates to a comprehensive verification of potential flooding signals identified from SAR backscatter time series.

[0008] This invention discloses a system capable of mapping rice planting intensity and planting and harvesting dates at a resolution of 10–20 m in any region of the world. The system relies solely on multi-source remote sensing data, meaning it eliminates the need for ground-based data such as field reference samples and statistical data. Therefore, it saves labor resources required for field sample collection, especially on large spatial scales. Specifically, it extracts potential planting signals from paddy fields using active microwave remote sensing data, and verifies the authenticity of these signals and estimates harvesting dates using optical, infrared, and passive microwave remote sensing data.

[0009] This invention uses troughs in the Sentinel-1C band SAR backscatter coefficient time series to detect rice transplanting time. Troughs in the SAR backscatter time series can be caused not only by rice transplanting but also by post-rainstorm flooding, vegetation cover changes, soil moisture or temperature variations, and unremoved noise related to changes in incident angle and topography. To compensate for this, this invention performs a six-dimensional test on all troughs. These dimensions are: 1) distance to adjacent troughs, 2) trough prominence, 3) trough width, 4) trough sharpness, 5) trough location time, and 6) trough value. The effective location time of the troughs is constrained using enhanced vegetation index (EVI) climatology retrieved from optical remote sensing and land surface temperature (LST) retrieved from multi-source infrared remote sensing. By revealing the influence of intra-annual variations in surface soil moisture and LST on the threshold, the spatial distribution pattern of the threshold trough values ​​(i.e., the SAR backscatter threshold used to identify the true flooding state) is mapped.

[0010] Therefore, this invention is an improved SAR-based high-resolution rice mapping system that is versatile and automated. Furthermore, it does not rely on field reference samples or statistical data that are typically essential to existing systems. Attached Figure Description

[0011] This patent or application document contains at least one color drawing. A copy of this patent or application publication with color drawings will be provided upon request and payment of the necessary fees by the Patent Office.

[0012] The foregoing and other objects and advantages of the invention will become more apparent when considered in conjunction with the following detailed description and accompanying drawings, wherein like reference numerals denote like elements in the various views, and wherein:

[0013] Figure 1 This is a flowchart illustrating the operation of one or more processors implementing the system according to the present invention;

[0014] Figure 2 This is a four-year composite map of rice planting intensity in the Asian monsoon region (excluding Northwest China) drawn using this invention;

[0015] Figure 3 It is made according to the present invention and Figure 2 Maps of the same area show the ordinal numbers of the 6-day cycle each year when the first rice crop is planted; and

[0016] Figure 4 It is made according to the present invention and Figure 2 The map of the same area shows the ordinal number of the 6-day cycle each year during the first rice harvest. Detailed Implementation

[0017] This invention relates to a system executed by one or more processors for detecting rice transplanting dates based on time series of valley backscatter coefficients in synthetic aperture radar (e.g., Sentinel-1C-band SAR). Considering that valleys in the SAR backscatter time series can be caused not only by rice transplanting but also by post-rainstorm flooding, vegetation cover changes, soil moisture or temperature variations, and unremoved noise associated with changes in incident angle and topographic effects, the system of this invention performs a six-dimensional test on all valleys. These dimensions are: 1) distance to adjacent valleys, 2) valley prominence, 3) valley width, 4) valley sharpness, 5) valley location time, and 6) valley value. Valid valley location times should be determined using enhanced vegetation index (EVI) climatology retrieved from optical remote sensing and land surface temperature (LST) retrieved from multi-source infrared remote sensing. By revealing the impact of intra-annual variations in surface soil moisture and LST on water depth thresholds, the spatial distribution pattern of threshold valley values ​​(i.e., SAR backscatter thresholds used to identify flooding conditions) can be mapped.

[0018] Figure 1 The diagram shows a technical flowchart of the operation of one or more processors of the system according to the present invention, wherein detailed steps are discussed below in

[0019] -

[0051] (in Figure 1 The steps (marked S1-S7) are merely specific embodiments of the present invention and are not intended to limit the scope of protection of the present invention. Any modifications or substitutions that are obvious to those skilled in the art are considered to be within the scope of protection of the present invention as defined by the appended claims.

[0019] S1. Sentinel-1 SAR Data Processing

[0020] The European Space Agency's (ESA) Sentinel-1 Earth observation mission consists of two satellites: Sentinel-1A (April 2014 to present) and Sentinel-1B (April 2016 to January 2022). Both satellites have a 12-day revisit period, resulting in two observations within the same 12-day period. For most of the world, the C-band synthetic aperture radar (SAR) instruments on both satellites operate in interferometric wide stripe (IW) mode. For Sentinel-1, IW products are typically available in single co-polarization (VV, vertical transmit, vertical receive) and cross-polarization (VH, vertical transmit, horizontal receive) modes. Because vertical-horizontal (VH) backscattering is superior to VV backscattering in detecting flooding conditions in paddy fields (Xu et al., 2023; Yang et al., 2021), in step 10, the present invention uses the VH band in a calibrated, orthogonally corrected ground range detection (GRD) scenario. To mitigate the prevalent noise in SAR images, it is recommended to apply additional boundary noise removal, single-temporal Lee-Sigma speckle filtering (Lee et al., 2009), and incident angle and radiative topography normalization (Mullissa et al., 2021; Vollrath et al., 2020) in step 12.

[0021] Although the pixel spacing of the IW stripes is 10×10m, the actual spatial resolution is 20×22m, while the kernel size used in the speckle filtering is 7×7. Therefore, embodiments of the present invention aggregate the preprocessed 10m resolution VH backscatter images to 20m to reduce uncertainty. Furthermore, to further reduce uncertainty, it is recommended to aggregate the VH backscatter to a 12-day resolution by calculating the median of two observations. Then, in step 14, a median absolute deviation (MAD) filter (Leys et al., 2013; Xu et al., 2018) and a Savitzky-Golay filter with a window size of 5 periods (60 days) should be used to remove potential outliers in the VH backscatter time series. Additionally, the Savitzky-Golay filter should be set to order 2 (Pan et al., 2021; Shen et al., 2023) to smooth the time series.

[0022] S2 surveys potential farmland

[0023] Because it is difficult to distinguish between paddy fields, constructed wetlands (Chen et al., 2013), and natural wetlands (Han et al., 2021; Zhou et al., 2016) on a large scale, according to the present invention, in step 16, all potential farmland is first mapped by merging two recognized state-of-the-art farmland maps. The first is the Global Land Analysis & Discovery (GLAD) Global 30m Farmland Map Standard 2016–2019 (Potapov et al., 2022); the second is ESAWorldCover 10m (2021) v200 (Zanaga et al., 2022). Both datasets should be resampled from their original resolution (10 or 30m) to 20m to match the resolution of the processed backscattered data. Furthermore, to avoid omissions, it is recommended to apply 2-D sequential statistical filtering to 20m pixels that are not classified as potential farmland. If at least three potential farmland pixels exist within a 3×3 window around the pixel, they are reclassified as potential farmland.

[0024] S3 utilizes multi-source surface temperature to extract data on the rice growing season.

[0025] Daily minimum surface temperature (LST) dailymin This is commonly used to determine the growing season of rice (Zhang et al., 2017). Firstly, because the LST from morning to 1:30 AM is approximately the LST. dailymin Therefore, in step 18, the system of this invention uses the MYD11A1.061Aqua nighttime LST dataset (Wan et al., 2021) to calculate the 12-day average LST. LST based on MODIS dailymin The time series should then undergo outlier removal and time filtering. This is because some cloud-prone areas lack LST (Long-Term Stamping). dailymin Due to the data period, this invention uses global 5km / h all-sky LST data (GLASS GHA-LST) (Jia et al., 2023). The LST data in this dataset... dailymin It can be used to generate LSTs with a 6-day resolution, such as those averaged over the period of interest. dailymin Climatology. The next step in step 20 is to apply the GLASS-based LST... dailymin Climatology resolution was reduced from 5km to 1km. This will be based on MODIS's 1km 1st... dailymin Time series data are aggregated to a 5km resolution, and then LST data are extracted from each 1km grid. dailymin Time series and corresponding LST in a 5km grid dailymin Time series regression was performed. The regression coefficients derived from each 1km grid were then used for GLASS-based LST. dailymin Downscaling in climatology.

[0026] Although previous studies extracted the start and end dates of the rice's hot growing season as the LST (longest temperature range) above 5°C in four consecutive 6-day intervals. dailymin The first and last dates (Zhang et al., 2017) are used, but considering the potential uncertainty of LST retrieved from remote sensing, the system of the present invention will use the LST in step 22. dailymin The threshold is adjusted to 0℃. Under no circumstances should rice seedlings be transplanted into water that may freeze overnight. By extracting the first and last 6-day cycles without freezing, the number of days between the first and last 12-day cycles without freezing can be determined as the annual hot growing season for rice.

[0027] S4 uses Sentinel-2-based EVI to extract rice heading stage.

[0028] The Sentinel-2 mission also includes two satellites, Sentinel-2A (June 2015 to present) and Sentinel-1B (March 2017 to present). These satellites carry a multispectral instrument (MSI). After processing poorly retrieved data (e.g., clouds and shadows) from the coordinated Sentinel-2 Level-2A data on Google Earth Engine (GEE), the present invention calculates the maximum EVI (Effective Virtualization Index) during the period of interest at step 24. max EVI should be excluded. max Potential farmland pixels below 0.4 (Han et al., 2021; Kontgis et al., 2015).

[0029] Next, for the remaining potential farmland pixels, the EVI climatology for the period of interest can be generated in step 26 by calculating the multi-year maximum EVI during each 10-day period. Here, the temporal resolution of the EVI climatology is set to 10 days because the revisit time for each individual Sentinel-2 satellite is 10 days, and the revisit time for the combined constellation is 5 days. This invention does not recommend taking into account the frequent cloud contamination in some tropical / subtropical regions when calculating the seasonal variation of the EVI for each individual year. The assumption for generating the EVI climatology is that the crop calendar rarely changes significantly over short time periods, such as less than 5 years.

[0030] The multi-year EVI climate can then be used in step 28 to extract the potential heading date (i.e., the heading date of rice, if rice is present) for each potential farmland pixel. Here, it is assumed that the rice heading date is close to a day of the year when the EVI peaks (DOY) (Zheng et al., 2016). Because the derived EVI climate typically contains noise or missing values ​​due to interannual variability in pollution and vegetation cover, but is highly periodic, this invention proposes the application of HANTS (Harmonic Analysis of Time Series) filtering, which is designed for vegetation index smoothing and phenological extraction (Roerink et al., 2000). In this invention, at most three planting seasons per year (i.e., three plantings) are considered, and the “number of frequencies considered above zero (nf)” is set to 3. Because the climatology is generated from the multi-year maximum EVI value, high and low outliers should not be excluded. The "fitting error tolerance" (fet) can be set to 5, the "degree of over-determinedness (dod)" can be set to 1, and the "suppression of small positive numbers with high amplitude (delta)" can be set to 0.1. To eliminate bad filtering results, HANTS filtering should not be applied to pixels with more than 18 periods (half of all 36 6-day periods) that lack EVI data. Furthermore, after HANTS filtering, the invention relates to calculating the Pearson correlation coefficient (r) between the original time series and the filtered time series, and if r is less than 0.6, the filter result should be considered invalid. Then, for pixels with valid filter results, the latter half of the filtered time series should be added to the beginning of the filtered time series, and the first half of the filtered time series should be added to its end to reconstruct a 10-day resolution EVI time series covering two full years.

[0031] This invention relates to peak search in the reconstructed EVI time series at step 30. The minimum distance between peaks should be set to 9 periods (≈90 days), the minimum width of a peak at a semi-prominent point should be set to 3 periods (≈30 days), and the minimum prominence of a peak should be set to 3 / 4 times the standard deviation of the reconstructed EVI. After the peak search, only peaks within one year (positions between periods 19 and 54) are retained. Considering that HANTS filtering can slightly alter the position of EVI peaks, according to an embodiment of the invention, a simple algorithm is applied to calibrate the derived peak positions. For each peak, the 10-day period with the highest EVI value is selected by comparing the original EVI values ​​from 2 periods (≈20 days) before the peak position to 2 periods after the peak position. If this highest EVI value is higher than both the median and mean of the original EVI time series and is also greater than 0.4, then the corresponding 10-day period can be considered the calibrated peak position. Specifically, if the distance between any two peaks after calibration is less than 9 periods (≈90 days), calibration should not be applied. Finally, the calibrated peak position should be adjusted to fit within the range between the first and 36th 10-day periods. For the few pixels lacking robust retrieval of peak position due to missing EVI data for more than 18 periods or invalid HANTS filtering, the potential heading date is estimated by sequential search and average neighbor robust retrieval (Chen et al., 2019). The total planting intensity in potential farmland pixels can be estimated based on the number of EVI peaks.

[0032] To match the temporal resolution of the processed SAR backscatter, the potential heading date at a 10-day resolution needs to be converted to DOY, and then to the ordinal number of the 12-day period within the period of interest. In step 32, taking into account the minimum duration of the vegetative growth stage and reproductive-to-maturity growth stage of rice in temperate or subtropical regions, and the basic cumulative temperature for rice growth (Liu et al., 2013; Zhang and Xu, 2012; Zhang et al., 2014), potential heading dates not exceeding 6 periods (≈72 days) after the start of the growing season and potential heading dates not exceeding 2 periods (≈24 days) before the end of the growing season are assumed not to belong to or indicate rice.

[0033] S5 extracts potential rice transplanting time by examining SAR backscatter valleys.

[0034] To extract potential rice transplanting signals, step 34 of this invention involves searching for valleys within the processed SAR backscatter time series. The minimum peak width of the semi-prominent point should be set to 2 periods (≈24 days) because the flooding period for rice, especially irrigated rice, is typically much longer than 24 days (Hashemi et al., 2022; Wei et al., 2023; Xu et al., 2023), while the low backscatter period caused by dryland rainstorms or random noise can hardly exceed 24 days (Najibi and Devineni, 2018). Similar to peak search in the EVI time series, considering the noise in the SAR data, the significance of the VH backscatter valley is forced to be higher than 3 / 4 times the standard deviation of the entire time series, and also greater than 1 dB in step 34. Furthermore, in step 34, valleys with backscatter values ​​higher than the mean or median of the entire SAR backscatter time series should be excluded. Because VH backscattered valleys caused by rice transplanting are generally sharper than those caused by changes in soil moisture, soil temperature, or random noise, one embodiment of the invention retains only peaks in step 34 whose sharpness is greater than 3 / 4 times the standard deviation of the entire time series and also greater than 1 dB, with the sharpness of the peaks calculated using Equation 1.

[0035]

[0036] In Equation 1, BS represents the SAR backscatter time series, i represents the ordinal number of a valley among all valleys, and Loc... i This indicates the position of the i-th valley.

[0037] The location and timing of backscattered valleys are also important. First, valleys located during the off-season should be removed in step 36. This step eliminates most backscattered valleys caused by very low vegetation cover, depleted soil moisture, and snow or ice cover. Second, for each potential heading stage or DOY (Dovey Year)... EVlmax In step 38, it is assumed that the corresponding potential transplanting period is in DOY. EVlmax -120 and DOY EVlmax Between -24. This time window should be further constrained by the location time of the previous potential heading date. For pixels with single-season or double-season planting, the lower limit of the time window should be higher than the previous DOY. EVlmax At least two cycles later (≈24 days). Furthermore, for pixels with three peaks in EVI climatology, because it is difficult to distinguish between triple cropping and three crops within two years, the lower limit of the time window is forced to be longer than the previous DOY. EVlmaxAt least 24 days later. Furthermore, because rice cannot overwinter like wheat, the lower limit of the allowed time window is lower than or equal to the start of the growing season. Third, in each time window used for transplanting signal detection, if more than one backscattering valley exists, embodiments of the invention retain only the sharpest one.

[0038] After extracting all effective valleys within the SAR backscatter time series, if the sharpness of the sharpest valley during, for example, the period from 2018 to 2021, is less than 3 / 4 of the standard deviation of the entire backscatter time series, or less than 1.8 dB (1 dB × (3 / 4)). 2 If all valleys are not considered potential rice planting signals, then in step 40, the VH backscatter values ​​(i.e., valley values) for all potential planting periods are output, which need to be further filtered to remove false ones (see the method in Section 6). Therefore, this invention relates to performing a six-dimensional test on all valleys in a SAR backscatter time series, which refers to testing the following: 1) distance to adjacent valleys; 2) valley prominence; 3) valley width; 4) valley sharpness; 5) valley location time; and 6) valley value.

[0039] S6 mapping threshold VH backscattering to identify true flooding status

[0040] Although it is generally accepted that the VH backscattering values ​​in submerged land are lower than those in unsubmerged land, and that the threshold VH backscattering value used to identify submerged states varies across different regions (Han et al., 2021), the spatial pattern of threshold VH backscattering has not yet been revealed. To fill this gap, this invention relates to mapping threshold Sentinel-1 VH backscattering at a resolution of 1 / 12°.

[0041] First, all VH backscatter values ​​during the potential transplanting period are collected within each 1 / 12° resolution grid. This is because there are 450... 2=202,500 pixels with a resolution of 20 meters (1 / 5400°), so the amount of data in each set can typically support unsupervised classification based on a Gaussian Mixture Model (GMM) (step 42). Here, it is assumed that the VH backscatter of the actual rice transplanting period is generally low and follows a Gaussian distribution, while those significant valleys in the VH backscatter time series that are unrelated to rice transplanting are generally high and follow another Gaussian distribution. Therefore, the GMM can simulate the mixed distribution of VH backscatter of the actual and spurious rice transplanting periods (Bishop and Nasrabadi, 2006). By estimating the parameters of the GMM using the Expectation-Maximization (EM) algorithm (step 44) (Moon, 1996), the threshold backscatter value can be estimated, and unsupervised classification between the two classes of VH backscatter valleys can be achieved. This method has previously been used to separate SAR backscatter in flooded farmland from SAR backscatter in unflooded land in small areas during the summer (Guan et al., 2023), but it has been modified in this invention and is applicable to large-scale and year-round flooding condition identification tasks.

[0042] In addition to adding the steps in S5, it is also meaningful to perform rigorous screening and spatial filtering on the 1 / 12° scale VH backscattering threshold estimation. First, the VH threshold should not be estimated for 1 / 12° grids with insufficient (<104) VH backscattering valley values ​​at the potential transplanting date. Second, according to previous studies (Bauer-Marschallinger et al., 2021; Xu et al., 2023), the fraction of VH backscattering in any group should not be less than 5%, and the mean of VH backscattering at the actual rice transplanting date should be in the range of [-28 dB, -16 dB]. Third, separation cannot be considered successful if the peak of a mixture component is located inside another component (Guan et al., 2023). Fourth, after calculating precision and recall according to Equation 2, models with biased predictions (precision-recall > 0.05) or low precision (precision+recall < 1.6) should be discarded; Fifth, when calculating the VH backscatter score classified as a true rice planting signal in each 1 / 12° grid, according to the present invention, VH backscatter threshold estimates that produce abnormally high / low scores are removed.

[0043] Specifically, in step 46, after applying a 2-D moving median filter with a window size of 18×18 grids (1.5°×1.5°) to the fractional plot of the actual transplanting signal covering the moving mean (MA), VH thresholds in grids where the difference between the original fractional estimate and the spatially smoothed fraction exceeds 0.2 should be removed. After applying a 2-D moving median filter with a window size of 1.5°×1.5° to all valid VH threshold estimates, the present invention relates to interpolating the smoothed result using biharmonic spline interpolation in step 48.

[0044] Precision=User's accuracy(UA)=1-commission error(CE)=TP / ((TP+FP));

[0045] Recall = Producer's accuracy (PA) = 1 - omission error (OE) = TP / ((TP+FN)); Equation 2

[0046] In Equation 2, TP, FP, and FN represent the number of true positives, false positives, and false negatives, respectively.

[0047] Next, the spatially interpolated threshold VH backscatter map can be calibrated by exploring the potential driving mechanisms in the spatial variation of threshold VH backscatter. It is known that, in addition to inundation, variations in surface soil moisture (SSM) can also lead to variations in VH backscatter. Low SSM is generally associated with low SAR backscatter (Hoskera et al., 2020; Sekertekin et al., 2020). Therefore, if a very low and sharp SSM valley exists during the potential transplanting period, it will produce a VH backscatter valley similar to the inundation signal. In this invention, an L-band passive microwave-based SSM product, namely the SMAPL3Radiometer Global Daily 9km Soil Moisture (SPL3SMP_E.005) (O'Neill et al., 2023), was also used. After performing outlier removal and Savitzky-Golay filtering on the 1 / 12° resolution SSM time series during 2018–2021 at step 50, key features (values, sharpness, and saliency) of the SSM valley can be calculated at step 52 for the potential planting period of all sub-pixels (20m resolution pixels). In addition to these features, the climatic background (P) should also be included. mean P cv LST mean LST std P is used as a predictor of the VH backscattering threshold. Here, P mean and P cv The mean annual precipitation and coefficient of variation for monthly precipitation climatology are obtained from the GPM IMERG final precipitation L3 half-hour 0.1° × 0.1° V06 dataset (Huffman et al., 2019), while LST mean LST stdThis represents the standard deviation of the average annual LST and monthly LST climatology. By training ten-fold random forest (RF) models in step 54 and applying these models in step 56, errors in the threshold VH backscatter map caused by interpolation in areas with insufficient data can be reduced.

[0048] S7 Rice Planting Intensity and Calendar Mapping

[0049] Using a 0.1° threshold VH backscatter map to identify the true flooding state, false planting periods can be removed from all potential planting periods extracted in S5. In step 58, the remaining true planting periods can then be used to plot the annual rice planting intensity. In an exemplary embodiment of the invention, a 2-D mode filter with a window size of 3×3 is applied to the annual rice planting intensity map to further reduce random errors (step 60). The invention also relates to developing a composite map every four years in step 62, which can be generated by first calculating the median of the rice planting intensity over the four years and then rounding the median (e.g., a composite of 0, 0, 1, 1 is 1; a composite of 0, 0, 0, 1 is "indeterminate"; and a composite of 0, 0, 0, 0 is 0). Specifically, in pixels where the rice planting intensity is lower than the total planting intensity estimated in S4, intercropping of rice with other crops is anticipated.

[0050] For pixels where rice is planted, the median DOY of all transplanting dates over a 4-year period is extracted for each potential rice growing season derived from S4. Then, in step 64, the heading date of that growing season is extracted (see S4), where the heading date is the date on which 50% of the plants produce seed heads in a typical year. Because the temporal resolution of the transplanting period extracted at this point is 12 days, the median DOY of the transplanting period can then be converted to an ordinal number for a 6-day cycle in a year.

[0051] Previous studies have shown that the time interval between harvest date and major rice-producing countries is 54 ± 10 days (More et al., 2016). Furthermore, the standard deviation of the time interval between harvest date and heading date is generally much smaller than the standard deviation of the time interval between harvest date and transplanting date (More et al., 2016; Zhao et al., 2023a). Therefore, the harvest date for rice in each growing season can be estimated to be 54 days later than the heading date. In all the identified paddy fields, some were directly sown rather than transplanted. For these fields, the rice planting date can also be considered to be approximately at the location time of the VH backscatter valley, as this is when the surface water content exceeds the crop vegetation (Hashemi et al., 2022).

[0052] 8. Application Examples of the Invention

[0053] As an example, this invention was applied to mapping rice planting intensity and calendar data during the Asian monsoon season from 2018 to 2021. A four-year composite map of rice planting intensity is shown below. Figure 2 As shown, this map is of monsoon Asia, where planting intensity is geographically located using color. Specifically, at the pixel level, areas without paddy fields or with uncertain locations are clearly shown; single-season paddy fields are shown in green, double-season paddy fields in orange, and triple-season paddy fields in red. The transplanting and harvesting time for the first season of rice (in 6 days) is shown in... Figure 3 and Figure 4 In the middle, they have the same basic color labels, but the difference is that they have more refined levels.

[0054] In summary, this invention is an improvement on existing SAR-based rice mapping because it is versatile, automated, and high-resolution. Therefore, it performs well in any region of the world, and it does not require field reference samples or statistical data, which are typically essential for existing systems.

[0055] References

[0056] References cited in this application are incorporated herein by reference in their entirety, as follows:

[0057] Amalia, RR, & Kadir, K. (2019). Improving Paddy Statistics in Indonesia: Developing Crop Cutting Survey Using Area Sampling Frame. In: University Library of Munich. Germany

[0058] Bauer-Marschallinger, B., Cao, S., Navacchi, C., Freeman, V., ReuB, F., Geudtner, D., Rommen, B., Vega, FC, Snoeij, P., Attema, E., Reimer, C., & Wagner, W. (2021). The normalized Sentinel-1 Global Backscatter Model, mapping Earth'sland surface with C-band microwaves. Scientific Data, S, 277

[0059] Bishop,C.M.,&Nasrabadi,N.M.(2006).Pattern recognition and machinelearning.Springer

[0060] Carlson,K.M.,Gerber,J.S.,Mueller,N.D.,Herrero,M.,MacDonald,G.K.,Brauman,K.A.,Havlik,P.,O’Connell,C.S.,Johnson,J.A.,Saatchi,S.,&West,P.C.(2017).Greenhouse gas emissions intensity of global croplands.Nature ClimateChange,7,63-68

[0061] Chen,H.,Zhu,Q.a.,Peng,C.,Wu,N.,Wang,Y.,Fang,X.,Jiang,H.,Xiang,W.,Chang,J.,Deng,X.,&Yu,G.(2013).Methane emissions from rice paddies naturalwetlands,lakes in China:synthesis new estimate.Global Change Biology,19.19-32

[0062] Chen,Y.,Feng,X.,Fu,B.,Shi,W.,Yin,L,&Lv,Y.(2019).Recent GlobalCropland Water Consumption Constrained by Observations.Water ResourcesResearch 55,3708-3738

[0063] Ding,Y.,Wang,W.,Zhuang,Q.,&Luo,Y.(2020).Adaptation of paddy rice inChina to climate change:The effects of shifting sowing date on yield andirrigation water requirement.Agricultural Water Management y 228,105890

[0064] Dong,J.,&Xiao,X.(2016).Evolution of regional to global paddy ricemapping methods:A review.ISPRS Journal of Photogrammetry and Remote Sensing,119,214-227

[0065] Guan,H.,Huang,J.,Li,L.,Li,X.,Miao,S.,Su,W.,Ma,Y.,Niu,Q.,&Huang,H.(2023).Improved Gaussian mixture model to map the flooded crops of VV and VHpolarization data.Remote Sensing of Environment y 295,113714

[0066] Han,J.,Zhang,Z.,Luo,Y.,Cao,J.,Zhang,L.,Cheng,F.,Zhuang,H.,Zhang,J.,&Tao,F.(2021).NESEA-RicelO:high-resolution annual paddy rice maps forNortheast and Southeast Asia from 2017to 2019.Earth Syst.Set Data,13,5969-5986

[0067] Han,J.,Zhang,Z.,Luo,Y.,Cao,J.,Zhang,L.,Zhuang,H.,Cheng,F.,Zhang,J.,&Tao,F.(2022).Annual paddy rice planting area and cropping intensity datasetsand their dynamics in the Asian monsoon region from 2000to2020.AgriculturalSystems,200,103437

[0068] Hashemi,M.G.Z.,Abhishek,A.,Jalilvand,E.,Jayasinghe,S.,Andreadis,K.M.,Siqueira,P.,&Das,N.N.(2022).Assessing the impact of Sentinel-1 derivedplanting dates on rice crop yield modeling.International Journal of AppliedEarth Observation and Geoinformation y 114,103047

[0069] Hoskera,A.K.,Nico,G.,Irshad Ahmed,M.,&Whitbread,A.(2020).Accuraciesof Soil Moisture Estimations Using a Semi-Empirical Model over Bare SoilAgricultural Croplands from Sentinel-l SAR Data.In,Remote Sensing

[0070] Huffman,G.J.,Stocker,E.F.,Bolvin,D.T.,Nelkin,E.J.,&Tan,J.(2019).GPMIMERG Final Precipitation L3 Half Hourly 0.1degree x 0.1degree V06.InG.E.S.D.a.I.S.C.G.DISC)(Ed.).Greenbelt,MD,USA:Goddard Earth Sciences Data andInformation Services Center (GES DISC)

[0071] Jia,A.,Liang,S.,Wang,D.,Ma,L.,Wang,Z.,&Xu,S.(2023).Global hourly,5 km,all-sky land surface temperature data from 2011 to 2021 based onintegrating geostationary and polar-orbiting satellite data.Earth Syst.SetData,75,869-895

[0072] Kontgis,C.,Schneider,A.,&Ozdogan,M.(2015).Mapping rice paddy extentand intensification in the Vietnamese Mekong River Delta with dense timestacks of Landsat data.Remote Sensing of Environment,169,255-269

[0073] Laborte,A.G.,Gutierrez,M.A.,Balanza,J.G.,Saito,K.,Zwart,S.J.,Boschetti,M.,Murty,M.V.R.,Villano,L.,Aunario,J.K.,Reinke,R.,Koo,J.,Hijmans,R.J.,&Nelson,A.(2017).RiceAtlas,a spatial database of global rice calendarsand production.Scientiffc Data,4,170074

[0074] Leys,C.,Ley,C.,Klein,0.,Bernard,P.,&Licata,L.(2013).Detectingoutliers:Do not use standard deviation around the mean.use absolute deviationaround the median.Journal of Experimental Social Psychology,49.764-766

[0075] Li,H.,Wang,X.,Wang,S.,Liu,Y.,Liu,Z.,Chen,S.,Wang,Q.,Zhu,T.,Wang,L.,&Wang,L.(2023).ChinaRiceCalendar-Seasonal Crop Calendars for Early,Middle,andLate Rice in China.Earth Syst.Sci.Data Discuss. y 2023,1-22

[0076] Liu,L.,Wang,E.,Zhu,Y.,Tang,L.,&Cao,W.(2013).Effects of warming andautonomous breeding on the phenological development and grain yield ofdouble-rice systems in China.Agriculture y Ecosystems&Environment y 165,28-38

[0077] Liu,Z.,Hu,Q.,Tan,J.,&Zou,J.(2019).Regional scale mapping offractional rice cropping change using a phenology-based temporal mixtureanalysis.International Journal of Remote Sensing,40,2703-2716

[0078] Mishra,B.,Busetto,L.,Boschetti,M.,Laborte,A.,&Nelson,A.(2021).RICA:Arice crop calendar for Asia based on MODIS multi year data.InternationalJournal of Applied Earth Observation and Geoinformation y 103,102471

[0079] Moon,T.K.(1996).The expectation-maximization algorithm.IEEE Signalprocessing magazine,13,47-60

[0080] More,R.S.,Manjunath,K.R.,Jain,N.K.,Panigrahy,S.,&Parihar,J.S.(2016).Derivation of rice crop calendar and evaluation of crop phenometrics andlatitudinal relationship for major south and south-east Asian countries:Aremote sensing approach.Computers and Electronics in Agriculture,127,336-350

[0081] Murali Krishna,G.,Andrew,N.,Prasad,S.T.,&Amrendra,N.S.(2011).Mappingrice areas of South Asia using MODIS multitemporal data.Journal of AppliedRemote Sensing,5,053547

[0082] Najibi,N.,&Devineni,N.(2018).Recent trends in the frequency andduration of global floods.Earth Syst.Dynam. y 9,757-783

[0083] O′Neill,P.E.,Chan,S.,Njoku,E.G.,Jackson,T.,Bindlish,R.,Chaubell,J.,&Colliander,A.(2023).SMAP Enhanced L3 Radiometer Global and Polar Grid Daily9km EASE-Grid Soil Moisture,Version 5.In N.N.S.a.I.D.Center(Ed.).Boulder,Colorado USA:Active Archive Center

[0084] Pan,B.,Zheng,Y.,Shen,R.,Ye,T.,Zhao,W.,Dong,J.,Ma,H.,&Yuan,W.(2021).High Resolution Distribution Dataset of Double-Season Paddy Rice inChina.In,Remote Sensing

[0085] Potapov,P.,Turubanova,S.,Hansen,M.C.,Tyukavina,A.,Zalles,V.,Khan,A.,Song,X.-P.,Pickens,A.,Shen,Q.,&Cortez,J.(2022).Global maps of cropland extentand change show accelerated cropland expansion in the twenty-firstcentury.Nature Food,3,19-28

[0086] Roerink,G.J.,Menenti,M.,&Verhoef,W.(2000).Reconstructing cloudfreeNDVI composites using Fourier analysis of time series.International Journalof Remote Sensing,21,1911-1917

[0087] Sekertekin,A.,Marangoz,A.M.,&Abdikan,S.(2020).ALOS-2and Sentinel-1SAR data sensitivity analysis to surface soil moisture over bare andvegetated agricultural fields.Computers and Electronics in Agriculture,171,105303

[0088] Shen,R.,Pan,B.,Peng,Q.,Dong,J.,Chen,X.,Zhang,X.,Ye,T.,Huang,J.,&Yuan,W.(2023).High-resolution distribution maps of single-season rice in Chinafrom 2017 to 2022.Earth Syst.ScL Data,15,3203-3222

[0089] Singha,M.,Dong,J.,Zhang,G.,&Xiao,X.(2019).High resolution paddy ricemaps in cloud-proue Bangladesh and Northeast India using Sentinel-ldata.Scientific Data,(5,26

[0090] Song,P.,Mansaray,L.R.,Huang,J.,&Huang,W.(2018).Mapping paddy riceagriculture over China using AMSR-E time series data.ISPRS Journal ofPhotogrammetry and Remote Sensing,144,469-482

[0091] Sun,C.,Zhang,H.,Xu,L.,Ge,J.,Jiang,J.,Zuo,L.,&Wang,C.(2023).Twenty-meter annual paddy rice area map for mainland Southeast Asia using Sentinel-1synthetic-aperture-radar data.Earth Syst.Set Data,75,1501-1520

[0092] Tian,G.,Li,H.,Jiang,Q.,Qiao,B.,Li,N.,Guo,Z.,Zhao,J.,&Yang,H.(2023).AnAutomatic Method for Rice Mapping Based on Phenological Features withSentinel-1 Time-Series Images.In,Remote Sensing

[0093] Wan,Z.,Hook,S.,&Hulley,G.(2021).MODIS / Aqua Land Surface Temperature / Emissivity Daily L3 Global 1km SIN Grid V061.In N.E.LP.D.A.A.Center(Ed.):NASAEOSDIS Land Processes Distributed Active Archive Center

[0094] Wang,X.,Folberth,C.,Skalsky,R.,Wang,S.,Chen,B.,Liu,Y.,Chen,J.,&Balkovic,J.(2022).Crop calendar optimization for climate change adaptation inrice-based multiple cropping systems of India and Bangladesh.Agricultural andForest Meteorology,315,108830

[0095] Wang,X.,Zhou,M.,Li,T.,Ke,Y.,&Zhu,B.(2017).Landuse change effects onecosystem carbon budget in the Sichuan Basin of Southwest China:Conversion ofcropland to forest ecosystem.Science of The Total Environment,609,556-562

[0096] Wei,J.,Cui,Y.,Luo,W.,&Luo,Y.(2022).Mapping Paddy Rice Distributionand Cropping Intensity in China from 2014to 2019with Landsat Images,EffectiveFlood Signals,and Google Earth Engine.In,Remote Sensing

[0097] Wei,J.,Cui,Y.,&Luo,Y.(2023).Rice growth period detection and paddyfield evapotranspiration estimation based on an improved SEBAL model:Considering the applicable conditions of the advectionequation.AgriculturalWater Management y 278,108141

[0098] Xu,S.,Zhu,X.,Chen,J.,Zhu,X.,Duan,M.,Qiu,B.,Wan,L,Tan,X.,Xu,Y.N.,&Cao,R.(2023).A robust index to extract paddy fields in cloudy regions from SARtime series.Remote Sensing of Environment,285,]A337 4

[0099] Xu,X.-Y.,Birol,F.,&Cazenave,A.(2018).Evaluation of Coastal Sea LevelOffshore Hong Kong from Jason-2 Altimetry.In,Remote Sensing

[0100] Yang,H.,Pan,B.,Li,N.,Wang,W.,Zhang,J.,&Zhang,X.(2021).A systematicmethod for spatio-temporal phenology estimation of paddy rice using timeseries Sentinel-1 images.Remote Sensing of Environment,259,1Y2394

[0101] Zanaga,D.,R.,V.D.K.,Daems,D.,De Keersmaecker,W.,Brockmann,C.,Kirches,G.,Wevers,J.,Cartus,0.,Santoro,M.,Fritz,S.,Lesiv,M.,Herold,M.,Tsendbazar,N.-E.,Xu,P.,Ramoino,F.,&Arino,O.(2022).ESA WorldCover 10m 2021v200.In E.S.Agency(Ed.)

[0102] Zhang,G.,Xiao,X.,Biradar,C.M.,Dong,J.,Qin,Y.,Menarguez,M.A.,Zhou,Y.,Zhang,Y.,Jin,C.,Wang,J.,Doughty,R.B.,Ding,M.,&Moore,B.(2017).Spatiotemporalpatterns of paddy rice croplands in China and India from2000to 2015.Science of The Total Environment,579,82-92

[0103] Zhang,G.,Xiao,X.,Dong,J.,Xin,F.,Zhang,Y.,Qin,Y.,Doughty,R.B.,&Moore,B.(2020).Fingerprint of rice paddies in spatial-temporal dynamics ofatmospheric methane concentration in monsoon Asia.Nature Communications,11,554

[0104] Zhang,J.,Wu,H.,Zhang,Z.,Zhang,L,Luo,Y.,Han,J.,&Tao,F.(2022).AsianRice Calendar Dynamics Detected by Remote Sensing and Their ClimateDrivers.In,Remote Sensing

[0105] Zhang,J.,&Xu,Y.(2012).Detecting Major Phenological Stages of RiceUsing MODIS-EVI Data and Symletll Wavelet in Northeast China.In,2012 2ndInternational Conference on Remote Sensing,Environment and TransportationEngineering(pp.1-4)

[0106] Zhang,S.,Tao,F.,&Zhang,Z.(2014).Rice reproductive growth durationincreased despite of negative impacts of climate warming across China during1981-2009.European Journal of Agronomy,54,70-83

[0107] Zhang,X.,Shen,R.,Zhu,X.,Pan,B.,Fu,Y.,Zheng,Y.,Chen,X.,Peng,Q.,&Yuan,W.(2023).Sample-free automated mapping of double-season rice in China usingSentinel-1 SAR imagery.Frontiers in Environmental Science,11

[0108] Zhao,X.,Nishina,K.,Akitsu,T.K.,Jiang,L,Masutomi,Y.,&Nasahara,K.N.(2023a).Feature-based algorithm for large-scale rice phenology detectionbased on satellite images.Agricultural and Forest Meteorology,329,109283

[0109] Zhao,X.,Nishina,K.,Izumisawa,H.,Masutomi,Y.,Osako,S.,&Yamamoto,S.(2023b).Monsoon Asia Rice Calendar:a gridded rice calendar in monsoon Asiabased on Sentinel-1 and Sentinel-2 images.Earth Syst.Set Data Discuss. y 2023,1-33

[0110] Zheng,H.,Cheng,T.,Yao,X.,Deng,X.,Tian,Y.,Cao,W.,&Zhu,Y.(2016).Detection of rice phenology through time series analysis of ground-basedspectral index data.Field Crops Research,198,131-139

[0111] Zhou, Y., Xiao, X., Qin, Y., Dong, J., Zhang, G., Kou, W., Jin, C., Wang, J., & Li,

[0112] While the invention has been explained with reference to certain embodiments, it should be understood that various modifications will become apparent to those skilled in the art upon reading the specification. Therefore, it should be understood that the invention disclosed herein is intended to cover such modifications that fall within the scope of the appended claims.

Claims

1. A system for general and automated high-resolution rice planting intensity and calendar mapping, the system comprising one or more processors for performing the following steps: Remotely acquire C-band synthetic aperture radar data to extract potential rice transplanting signals from paddy fields; Acquire optical, infrared, and passive microwave remote sensing data; The authenticity of the potential rice transplanting signal was verified using the remote sensing data. Planting and harvesting dates are estimated based on certified transplanting signals to determine calendar mapping; as well as The rice planting intensity and calendar mapping are generated based on the certified transplanting signals without using ground-based data such as field reference samples and statistics.

2. The system according to claim 1, wherein, Search for valleys in the time series of backscatter coefficients (backscattering) of satellite-based C-band synthetic aperture radar (SAR) to determine the rice transplanting period.

3. The system according to claim 1, wherein, The resolution of the survey is between 10 meters and 20 meters in any region of the world.

4. The system according to claim 2, wherein, The step of verifying the authenticity of the potential rice planting signal using the remote sensing data is used to remove potential flooding signals identified from the SAR backscatter time series.

5. The system according to claim 4, wherein, Valleys unrelated to rice in the SAR backscatter time series are removed by performing a six-dimensional test on all valleys. The six-dimensional test includes: 1) distance to adjacent valleys, 2) valley prominence, 3) valley width, 4) valley sharpness, 5) valley location and time, and 6) valley value.

6. The system according to claim 5, wherein, The effective SAR backscatter valley location time of rice paddies was determined using enhanced vegetation index (EVI) climatology derived from optical remote sensing data and land surface temperature (LST) retrieved from multi-source infrared remote sensing images.

7. The system according to claim 6, wherein, Spatial patterns of threshold valleys (i.e., threshold SAR backscatter values ​​used to identify true flooding conditions) were mapped using a Gaussian mixture model, and the influence of intra-annual variations in surface soil moisture and LST on this threshold was revealed.

8. The system according to claim 1, wherein, The satellite-based C-band synthetic aperture radar (SAR) backscattering coefficient data were acquired in the VH band from the European Space Agency's (ESA) Sentinel-1 Earth observation mission, which consists of two satellites: Sentinel-1A (April 2014 to present) and Sentinel-1B (April 2016 to January 2022).

9. The system according to claim 1, wherein, The noise prevalent in SAR images is mitigated by using single-temporal Lee-Sigma speckle filtering and additional boundary noise removal with incident angle and radiation topography normalization.

10. The system according to claim 1, wherein, The 10m resolution VH backscatter time series image was preprocessed by the following steps: (a) clustering to 20m, (b) clustering to 12-day resolution by calculating the median of two observations to reduce uncertainty, (c) then removing outliers from the VH backscatter time series using a median absolute deviation (MAD) filter and a Savitzky-Golay filter, and (d) smoothing the time series with a second-order filter.

11. The system according to claim 1, further comprising the following steps: (a) All potential farmland was mapped using a combination of the global 30m farmland map from Global Land Analysis and Discovery (GLAD) and the ESAWorldCover 10m map, and then (b) the two map datasets were resampled from their original resolution (10m or 30m) to 20m to match the resolution of the processed backscatter data.

12. The system according to claim 1, further comprising the following steps: Apply 2-D sequential statistical filtering to 20m pixels that are not classified as potential farmland, and reclassify them as potential farmland if there are at least three potential farmland pixels within a 3×3 window around the pixel.

13. The system according to claim 6, wherein, The 12-day average daily minimum surface temperature (LST) was calculated using the MYD11A1.061 Aqua nighttime LST dataset. dailymin Then, outlier removal and time filtering are performed, and global hourly 5km all-sky LST data (GLASS GHA-LST) is added.

14. The system of claim 13, wherein (a) the LST of the dataset dailymin Used to generate a 6-day resolution LST averaged over the period of interest. dailymin Climatological data, (b) the GLASS-based LST dailymin Climatology reduced from 5km to 1km resolution, (c) 1km LST based on MODIS dailymin The time series are aggregated to a 5km resolution, and then (d) LSTs in each 1km grid are... dailymin Time series relative to LST in the corresponding 5km grid dailymin The time series is regressed, and then (e) the derived regression coefficients in each 1km grid are used for GLASS-based LST. dailymin The shrinking of climatology.

15. The system according to claim 6, wherein, In the extraction of non-freezing (LST) dailymin The number of days between the first and last 6-day cycles when the temperature is above 0℃, and between the first and last 12-day cycles when the temperature does not freeze, is determined as the hot growing season for rice each year.

16. The system according to claim 6, wherein, The maximum EVI climate within the period of interest is generated by calculating the multi-year maximum EVI in every 10-day period using coordinated Sentinel-2 Level-2A data. Potential farmland pixels with EVImax below 0.4 are then excluded. For the remaining potential farmland pixels, a peak search algorithm is applied to the smoothed EVI climate to extract the potential heading date of rice in each potential farmland pixel.

17. The system according to claim 16, wherein, Except for pixels with more than 18 periods (half of all 36 6-day periods) without EVI data, Time Series Harmonic Analysis (HANTS) filtering was applied to smooth the derived EVI climate signal before it was used to extract the heading period.

18. The system according to claim 17, wherein, After the HANTS filtering, the Pearson correlation coefficient (r) between the original time series and the filtered time series is calculated, and if r is less than 0.6, the filtering result is considered invalid. as well as For pixels with effective filter results, the latter half of the filtered time series is added to the head of the filtered time series, and the first half of the filtered time series is added to its tail to reconstruct a 10-day resolution EVI time series covering two full years.

19. The system according to claim 16, wherein, During the peak search in the reconstructed EVI time series, the minimum distance between peaks was set to 9 periods (~90 days), the minimum width of peaks at the half-peak was set to 3 periods (~30 days), and the minimum peak prominence was set to 3 / 4 times the standard deviation of the reconstructed EVI; and After peak search, only peaks within one year (those located between the 19th and 54th cycles) are retained.

20. The system according to claim 19, wherein, An algorithm is applied to calibrate the estimated peak position, as shown below: For each peak, compare the original EVI values ​​from two periods (~20 days) before the peak position to two periods after the peak position; Select the 10-day period with the highest EVI value; If the highest EVI value is higher than the median and mean of the original EVI time series and is also greater than 0.4, the corresponding 10-day period is considered to be the calibrated peak position. If the distance between the two peaks after calibration is less than 9 periods (≈90 days), then calibration should not be applied. The calibrated peak position is adjusted to fit within the range between the 1st and 36th 10-day periods, where the potential heading date is estimated by sequential search and average neighbor robust search for pixels that have no robust retrieval of peak position due to missing EVI data for more than 18 periods or invalid HANTS filtering.

21. The system according to claim 16, wherein, To match the temporal resolution of the processed SAR backscatter, the potential heading period at a 10-day resolution is converted to a single day of the year (DOY), and then to the ordinal number of a 12-day period within the period of interest; and Among them, it is assumed that the potential heading period of no more than 6 cycles (≈72 days) after the start of the growing season and the potential heading period of no more than 2 cycles (≈24 days) before the end of the growing season do not belong to or are not indicative of rice.

22. The system according to claim 5, wherein, The criterion used to search for valleys within the processed SAR backscatter time series is the minimum peak width at the half-protrusion over two periods (≈24 days). The significance of the VH backscattering valley was more than 3 / 4 times the standard deviation of the entire time series and also greater than 1 dB, excluding valleys where the backscattering value was higher than the mean or median of the entire SAR backscattering time series; and Only those peaks whose sharpness, calculated using Equation 1, is greater than 3 / 4 times the standard deviation of the entire time series, and also greater than 1 dB, are retained, where Equation 1 is... Where BS represents the SAR backscatter time series, i represents the ordinal number of a valley among all valleys, and Loc i This indicates the position of the i-th valley.

23. The system according to claim 6, wherein, During peak search, valleys located outside the growing season are removed; For each potential heading stage or DOY EVlmax Due to the constraints of the previous potential heading period location and timing, it is assumed that the corresponding potential transplanting period is in DOY. EVlmax -120 and DOY EVlmax Between -24; For pixels with three peaks in EVI climatology, the lower time bound is forced to be at the second previous DOY. EVlmax At least 24 days afterwards; The lower limit of the allowable time window is lower than or equal to the start of the growing season; as well as In each time window used for rice transplanting signal detection, if there is more than one backscatter valley, only the sharpest one is retained.

24. The system according to claim 5, wherein, After extracting all valid valleys within the SAR backscatter time series, if the sharpest valley's sharpness is less than 3 / 4 of the standard deviation of the entire backscatter time series or less than 1.8 dB, then the valley should not be considered an indication of a potential rice transplanting signal; and Output the VH backscatter values ​​(i.e., valley values) for all potential transplanting periods.

25. The system according to claim 7, wherein, By (a) acquiring all VH backscatter values ​​during the potential transplanting period within each 1 / 12° resolution grid, (b) applying unsupervised classification based on a Gaussian mixture model (GMM), and (c) estimating the parameters of the GMM using the expectation-maximization (EM) algorithm, it is possible to estimate the threshold backscatter value and achieve unsupervised classification between two types of VH backscatter valleys, mapping the threshold Sentinel-1 VH backscatter time series at 1 / 12° resolution; and Specifically, for 1 / 12° grids with insufficient (<104) VH backscattering valley values ​​at the potential transplanting period, the VH threshold is not estimated; when it is below 5%, the VH backscattering fraction is not considered; and the mean VH backscattering at the actual rice transplanting period is in the range of [-28dB, -16dB]. Where the peak of a mixture component lies within another component, separation is not considered. After calculating precision and recall according to Equation 2, models with biased predictions ((|Precision-Recall|>0.05) or low precision (precision+recall<1.6) are discarded. Furthermore, when calculating the VH backscattering fraction classified as a true transplanting signal in each 1 / 12° grid, VH backscattering threshold estimates that produce abnormally high / low fractions are removed, where Equation 2 is... Precision=User's accuracy(UA)=1-commission error(CE)=TP / ((TP+FP)); Recall = Producer's accuracy (PA) = 1 - omission error (OE) = TP / ((TP+FN)) where TP, FP, and FN are the number of true positives, false positives, and false negatives, respectively.

26. The system according to claim 25, wherein, After applying a 2-D moving median filter 46 with a window size of 18×18 grids (1.5°×1.5°) to the graph of the fractions of the real transplanting signal covering the moving mean (MA), the VH threshold in the grids where the difference between the original fraction estimate and the spatially smoothed fraction 52 exceeds 0.2 is removed. as well as After applying a 2-D moving median filter with a window size of 1.5° × 1.5° to all valid VH threshold estimates, the results are smoothed using biharmonic spline interpolation.

27. The system according to claim 7, wherein, The threshold VH backscatter map of spatial interpolation is calibrated by exploring the potential driving mechanism in the spatial variation of threshold VH backscatter. L-band passive microwave-based SSM product, namely SMAP L3 global radiometer; Utilizing 9km of daily soil moisture (SPL3SMP_E.005); and After performing outlier removal and Savitzky-Golay filtering on the 1 / 12° resolution SSM time series, key features (values, sharpness, and prominence) of the SSM valleys were calculated for all potential planting periods at sub-pixels (20m resolution pixels).

28. The system according to claim 7, wherein, Climate background (Pmean, Pcv, LSTmean, LSTstd) are used as predictors of the VH backscattering threshold. Pmean and Pcv are the coefficients of variation of the mean annual and monthly precipitation climatology obtained from the GPM_IMERG final precipitation L3 half-hour 0.1° × 0.1° V06 dataset, while LSTmean and LSTstd are the standard deviations of the mean annual and monthly LST climatology. By training and applying ten-fold random forest (RF) models, errors in the threshold VH backscatter map caused by interpolation in regions with insufficient data were reduced.

29. The system according to claim 5, wherein, The true flooding state is identified by using the 0.1-degree backscattered map of the threshold VH, spurious information is removed from all extracted potential transplanting periods, and the remaining true transplanting periods are then used to plot the annual rice planting intensity.

30. The system according to claim 1, wherein, A 2-D mode filter with a window size of 3×3 was applied to the annual rice planting intensity map to further reduce random errors.

31. The system according to claim 1, wherein, Generate composite charts for each four-year period using the following method: First, calculate the median of rice planting intensity over the four-year period, and then round the median to the nearest whole number. For pixels where rice is planted, extract the median DOY of all transplanting dates over a 4-year period for each potential rice growing season, then extract the heading date for that growing season; and Since the time resolution of the rice transplanting period extracted at this time is 12 days, the median value DOY of the rice transplanting period is then converted into an ordinal number of a 6-day cycle in a year.

32. The system according to claim 1, wherein, The estimated harvest date for rice in each growing season is 54 days later than the heading date.

Images