A remote sensing method for extracting rice growing season based on temporal bidirectional detection

Abstract

This invention discloses a remote sensing method for extracting the rice growing season based on time-series bidirectional detection. The method includes constructing a remote sensing index time-series dataset, determining the effective planting season and planting system, extracting the Start of Season (SOS) of the growing season, identifying the rice season, and extracting key parameters of the rice growing season. This invention can accurately identify the rice season and obtain high-precision key phenological periods of rice even under complex planting systems. Compared with traditional remote sensing methods, this method can accurately identify the crop planting system in the study area, and the extraction accuracy of the rice season and its key phenological periods is higher. Furthermore, this method can achieve cross-scale transfer, has strong versatility, and is applicable to rice crops in other study areas. By accurately identifying pixels containing rice planting and their rice seasons, it can accurately monitor the growth and development process of rice crops and extract key phenological periods.

Description

TECHNICAL FIELD

[0001] The present application belongs to the technical field of remote sensing parameter inversion, and particularly relates to a rice growing season remote sensing extraction method based on time sequence bidirectional detection. BACKGROUND

[0002] Rice, as one of the three major crops in the world, is the main food source for more than half of the world's population. Its phenological phase directly affects resource utilization efficiency and grain yield, and is of great significance to agricultural precision management, food security and ecological system stability. Under the background of global climate change, accurate identification of rice phenology can effectively reduce the impact of extreme weather on yield and provide scientific support for agricultural adaptation strategies.

[0003] Traditional rice phenology acquisition relies on manual ground observation, which has the defects of high labor and time cost, non-uniform observation standards, and incomplete spatial coverage, and is difficult to meet the needs of large-scale monitoring. Remote sensing technology has become the mainstream method for crop phenology monitoring due to its macroscopic and timely advantages. By constructing a time series curve of vegetation index, combined with inflection point identification and threshold determination, the phenological parameters are extracted. In the prior art, the paper Luo Y, Zhang Z, Chen Y, et al. China Crop Phen1km: a high-resolution crop phenological dataset for three staple crops in China during 2000-2015 based on leaf area index (LAI) products [J]. Earth System Science Data, 2020, 12(1): 197-214. describes the use of LAI products combined with various filtering methods to extract crop phenology. The paper Xin Z, Kazuya N, Kawaguchi T A, et al. Feature-based algorithm for large-scale rice phenology detection based on satellite images [J]. Agricultural and Forest Meteorology, 2023, 329. combines optical and microwave remote sensing to improve spatial resolution, but these methods still have deficiencies in complex planting system areas. SUMMARY

[0004] The present application aims to provide a rice growing season remote sensing extraction method based on time sequence bidirectional detection, and solve the technical problems of planting season misjudgment, rice identification difficulty and low growth season parameter precision in the prior art for rice growing season extraction in a complex planting system region.

[0005] To solve the above technical problems, the present application adopts the following technical solutions:

[0006] A rice growing season remote sensing extraction method based on time sequence bidirectional detection, comprising the following steps:

[0007] Step S1: constructing a remote sensing index time sequence dataset, specifically comprising:

[0008] Step S1.1: extracting a cultivated land or farmland pixel set S in the research area based on a high spatial resolution land use classification product;

[0009] Step S1.2: after preprocessing the spectral reflectance remote sensing image covering the research area, calculating the original time sequence of normalized vegetation index NDVI, enhanced vegetation index EVI and land surface water index LSWI pixel by pixel, respectively denoted as NDVI original sequence S n , EVI original sequence S e , and LSWI original sequence S l ;

[0010] Step S1.3: using Spline interpolation method to fill in the missing values in S n , S e , and S l , to obtain the interpolated sequences S n ', S e ', and S l ';

[0011] Step S1.4: using Savitzky-Golay filtering algorithm to smooth the interpolated sequences S n ', S e ', and S l ', and the filtering process needs to include remote sensing index data of the previous December and the following January of the research year to reduce boundary effects, finally retaining the effective data of the research year to form a remote sensing index time sequence dataset D={D n , D e , D l}, wherein D n is the NDVI time sequence of the research year, D e is the EVI time sequence of the research year, and D l is the LSWI time sequence of the research year.

[0012] Step S2: determining the effective planting season and planting system, specifically comprising:

[0013] Step S2.1: Analyze the EVI time series D for the study year. e All vegetation index values ​​are sorted from smallest to largest, and the 80th percentile is taken as the peak threshold for the planting season. p ;

[0014] Step S2.2: Identify the EVI time series D for the research year e Curve and Th p The intersection points, from front to back, each pair of consecutive intersection points constitutes an initial planting season S. ini Let S be the initial set of planting seasons. ini-set ={S ini 1 ,S ini 2 ,...,S ini k}, where k is the initial number of planting seasons;

[0015] Step S2.3: For each initial planting season S ini Extract the corresponding EVI peak value, denoted as EVI peak value P. i And it is verified through the following three constraints: time constraint, amplitude constraint and interval constraint;

[0016] Step S2.4: Retain the initial planting season that satisfies all the above constraints, and denote it as the effective planting season S. val The effective planting season is S. val-set ={S val 1 ,S val 2 ,...,S val m}, where m is the number of effective planting seasons, which is the planting system of the pixel. m=1 is single-season planting, and m≥2 is multi-season planting.

[0017] Step S3: Extract the effective start date SOS of the planting season, specifically including:

[0018] Step S3.1: For each effective planting season S val Extract its time series D in the EVI of the study year e The maximum value of EVI in the data is denoted as the peak EVI for that planting season. max-v The corresponding date sequence is DOY. max-v ;

[0019] Step S3.2: For different planting systems, starting from the day sequence corresponding to the peak EVI of each effective planting season, traverse the EVI time sequence backward to determine the minimum EVI value. Then, calculate the extraction threshold based on this and detect forward. The first day sequence that reaches the extraction threshold is the start date SOS of that planting season.

[0020] Step S3.3: Denote the set of SOS for all valid planting seasons as SOS. -set ={SOS1,SOS2,...,SOSm}

[0021] Step S4: Identify the rice season: Using the SOS of each valid planting season as the time origin, define a 20-day time window before and after it. Count the number of days within the window where the LSWI is consistently higher than the NDVI or EVI of the same period. If the number of consecutive days is ≥8, then the planting season is determined to be the rice season.

[0022] Step S5: Extract key parameters of the rice growing season: For the rice season, based on the differences in planting systems, starting from the day sequence corresponding to the peak EVI of the rice season, traverse the EVI time sequence backward to determine the minimum EVI value; then calculate the tail extraction threshold based on this, and detect the end date EOS of the rice growing season backward. Calculate the growing season duration by the difference between EOS and SOS, thus completing the extraction of the rice growing season.

[0023] Further optimization is achieved by using the following formulas for calculating the vegetation index NDVI, enhanced vegetation index EVI, and land surface moisture index LSWI in step 1:

[0024] ;

[0025] ;

[0026] ;

[0027] in, and These are the near-infrared, red, blue, and short-wave infrared spectral bands, respectively.

[0028] Further optimization is achieved by defining the peak constraint condition in step S2.3 as follows:

[0029] Time constraint: Peak value P i The corresponding day sequence falls within the crop growth period when the average daily temperature is greater than 5℃ for 5 consecutive days. Tg Within. Crop growth period Tg The determination needs to be based on the daily temperature data of the study area, and the average temperature over five consecutive days needs to be calculated using the sliding window method. The point at which the average temperature first exceeds 5°C is recorded as [a certain value]. Tg The starting date is defined as the last time the average temperature exceeds 5°C. Tg The end date, forming a complete Tg Interval.

[0030] Amplitude constraint: Peak value P i Corresponding Enhanced Vegetation Index EVI Value > 0.4.

[0031] Interval constraint: If k≥2, the maximum EVI peak value within the year is taken as the main peak P. max The remaining peak values ​​are considered as secondary peaks P.suβ Then P suβ With P max The corresponding day sequence difference ΔDOY = |DOY suβ -DOY max 60 days.

[0032] Further optimization is achieved by defining the SOS signal extraction rules for different planting regimes in step S3.2 as follows:

[0033] Step S3.2.1: If it is a single-season planting, it includes:

[0034] Step S3.2.1.1: Using DOY max-v Starting from the beginning, traverse the EVI time series D of the study year backwards. e Find the minimum value and denote it as the minimum EVI for a single season. min-s ;

[0035] Step S3.2.1.2: According to formula Th s =EVI min-s +0.2×(EVI max-v -EVI min-s ), calculate the SOS extraction threshold Th s ;

[0036] Step S3.2.1.3: Using DOY max-v Using forward detection as a starting point, we can study the annual EVI time series D. e The first one equals Th s The date corresponding to the point is the start date of the growing season (SOS) for that effective planting season. s。

[0037] Step S3.2.2: If it is multi-season planting, it includes:

[0038] Step S3.2.2.1: For the first effective planting season S val 1 Extract its SOS1 according to the method in step S3.2.1;

[0039] Step S3.2.2.2: For the t-th effective planting season S in multi-season planting. val t If t≥2, then it includes:

[0040] Step S3.2.2.2.1: Find the (t-1)th effective planting season S val t-1 peak EVI max-v t-1 and the corresponding day sequence DOY max-v t-1 ;

[0041] Step S3.2.2.2.2: Traverse the EVI time series D of the study year e DOY max-v t-1 To DOYmax-v t Find the interval D within that interval. e The minimum value is denoted as the minimum EVI of the t-th effective planting season. min-m t ;

[0042] Step S3.2.2.2.3: According to formula Th st =EVI min-m t +0.2(EVI max-v t -EV min-m t ), calculate the SOS extraction threshold Th for the t-th planting season. st ;

[0043] Step S3.2.2.2.4: Using DOY max-v t Using forward detection as a starting point, we can study the annual EVI time series D. e The first one equals Th st The day sequence corresponding to the point is S. val t SOS at the start of the growing season t .

[0044] Further optimization, step SS4 specifically includes:

[0045] Step S4.1: Set the starting date of the growing season (SOS) for the i-th effective planting season. i i=1,2,...,m; Define "DOY=SOS" i -20” to “DOY=SOS” i +20” time window W i ;

[0046] Step S4.2: Rice determination window W in the i-th effective planting season i Inside, compare D day by day l With D n D l With D e The numerical values ​​of the values ​​satisfy D. l >D n Or D l >D e The number of consecutive days is denoted as the number of continuous flooding days D for the i-th window. i ;

[0047] Step S4.3: If D i If the value is ≥8, then SOS is determined. i Corresponding effective planting season S val i For rice season R i Conversely, it is a non-rice season; let the set of rice seasons be R. -set ={R1,R2,...,R n}, where n is the number of rice seasons, and n≤m.

[0048] For further optimization, in step S5, it specifically includes:

[0049] Step S5.1: For the j-th rice season R j , where j = 1, 2,..., n, extract its EVI peak value in the EVI time series D e of the research year, and denote it as the EVI peak value of the rice season EVI max-r j , and the corresponding day order is denoted as DOY max-r j。

[0050] Step S5.2: If R j is the last effective rice planting season within the year, that is, j = m, then it includes:

[0051] Step S5.2.1: Starting from DOY max-r j , traverse the EVI time series D e of the research year backward to find the minimum value, and denote it as the EVI minimum value at the end of the rice season EVI min-rear ;

[0052] Step S5.2.2: Calculate the EOS end extraction threshold Th end according to the formula Th min-rear = EVI max-r j + 0.2×(EVI min-rear - EVI end );

[0053] Step S5.2.3: Starting from DOY max- r j , detect the EVI time series D e of the research year backward. The day order corresponding to the first point equal to Th end is the end of the growing season EOS j of R j 1。

[0054] Step S5.3: If there are subsequent effective planting seasons after R j , that is, j < m, then it includes:

[0055] Step S5.3.1: Find the EVI peak value EVI j of the next effective planting season S val j+1 of R max-v j+1 and the corresponding day order DOY max-v j+1 ;

[0056] Step S5.3.2: Traverse the interval from DOY e to DOY max-r j in the EVI time series D max-v j+1 of the research year to find the minimum value of D e within this interval, and denote it as the EVI minimum value in the middle of the rice season EVI min-rmid ;

[0057] Step S5.3.3: According to formula Th end j =EVI min-rmid +0.2×(EVI max-r j -EVI min-rmid ), calculate the EOS tail extraction threshold Th for the j-th rice season. end j ;

[0058] Step S5.3.4: Using DOY max-r j Starting from the point of origin, detect D backwards. e The first one equals Th end j The day sequence corresponding to the point is R. j End of the growing season for EOS j。

[0059] Step S5.4: Unify the notation of R j The growing season ends for EOS j , j=1,2,...,n, calculate R according to the formula j The corresponding growing season duration T j =EOS j -SOS j。

[0060] Step S5.5: Integrate all R j SOS j EOS j T j This forms the key parameter set Para={Para1,Para2,...,Para} for the rice growing season. n}, where Para j ={SOS j EOS j ,T j}, complete the extraction during the rice growing season.

[0061] Further optimization involves step S3.2, which forward-detects the study year EVI time series D. e The specific method is as follows: If it is a single-season planting, from DOY max-v Initially, read D day by day in descending order. e The value, when the read value first equals Th s When the time is right, stop the detection and record the corresponding day sequence; if it is multi-season planting, start from DOY max-v t Initially, read D day by day in descending order. e The value, when the read value first equals Th st When the value is reached, stop the detection and record the corresponding day sequence.

[0062] Further optimization involves step S4.2, where the number of consecutive flooded days D... i The statistical rule is: in W i Inside, when D appears l >D n Or D l >D e The timer starts when the condition is met. If the condition is met on subsequent dates, the count continues. If the condition is not met on a certain date, the count is reset, and the final value is W. i The longest consecutive count of days is taken as D. i .

[0063] Further optimization involves backward detection of the study year EVI time series D in steps S5.2.3 and S5.3.4. e The specific method is: from DOY max-r j Initially, read D day by day in ascending order. e The value, when the read value first equals Th end or Th end j When the value is reached, stop the detection and record the corresponding day sequence.

[0064] Compared with the prior art, the beneficial effects of the present invention are as follows:

[0065] 1. More accurate identification of planting system: By using EVI quantile threshold and three constraints, "interference peaks" and non-growing season signals are effectively eliminated, and single-season / multi-season planting systems are accurately identified, solving the problems of misjudgment and omission of planting season in traditional methods.

[0066] 2. More reliable starting point extraction during the growing season: It adopts forward detection bidirectional logic based on the EVI peak value, and adjusts the calculation benchmark according to the differences in planting system to avoid SOS offset caused by unidirectional detection and improve the accuracy of starting point extraction.

[0067] 3. More efficient rice season differentiation: Utilizing the unique LSWI spectral response of rice during the flooding period of rice transplanting, the continuous feature determination within a 40-day time window solves the rice identification problem caused by "different spectra for the same object" and "same spectra for different objects".

[0068] 4. Enhanced applicability to large scales: Based on satellite remote sensing data and combined with MATLAB programs for automated processing, it can efficiently process large-scale data over long time periods and large areas, making it suitable for large-scale applications in regions with complex cropping systems. Attached Figure Description

[0069] Figure 1 This is a flowchart of the method described in this invention;

[0070] Figure 2This is a schematic diagram of the remote sensing index time series curve of rice planting pixels in Example 1; where (a) is the remote sensing index time series curve, (b) is a schematic diagram of determining the effective planting season, (c) is a schematic diagram of extracting the SOS of each effective planting season, (d) is a schematic diagram of determining whether each effective planting season is a rice season, (e) is a schematic diagram of extracting the key phenological period of each rice season, and (f) is a schematic diagram of determining the rice season using traditional remote sensing methods.

[0071] Figure 3 The diagram shows the time series curves of remote sensing indices for four types of planting systems; (a) is the time series curve of remote sensing indices for non-rice + rice, (b) is the time series curve of remote sensing indices for non-rice + non-rice, (c) is the time series curve of remote sensing indices for single-season rice, and (d) is the time series curve of remote sensing indices for single-season non-rice. Detailed Implementation

[0072] The invention will now be further described with reference to the accompanying drawings. The following embodiments are only used to more clearly illustrate the technical solution of the invention and should not be used to limit the scope of protection of the invention.

[0073] In this embodiment, Jiangsu Province was selected as the study area. Jiangsu is located in the central part of the coastal area of ​​China (30°45′~35°08′ north latitude, 116°21′~121°56′ east longitude), in the lower reaches of the Yangtze River and Huai River. The terrain is mainly plains, with most areas below 50 meters above sea level. It has a transitional climate from temperate to subtropical, with abundant rainfall and heat. The agricultural planting system is complex, including single-season rice, double-season rice, and rotation patterns such as winter wheat-rice and rapeseed-rice. It is an important rice-producing area in my country and is suitable as a verification area for this method.

[0074] like Figure 1 As shown, a remote sensing method for extracting rice growing season data based on temporal bidirectional detection includes the following steps:

[0075] Step S1: Construct a time series dataset of remote sensing indices.

[0076] Taking 2023 as an example, the required data was obtained by using the European Space Agency's land use classification product ESAWorldCover to obtain the set of arable land or farmland pixels S in Jiangsu Province.

[0077] Then, the spectral reflectance remote sensing image data of Sentinel-2 satellite (with a time resolution of 5 days) was preprocessed using the Google Earth Engine (GEE) cloud computing platform to obtain 42 composite median images for 10 days from December 1, 2022 to January 31, 2024, with a spatial resolution of 10 meters.

[0078] Next, based on the near-infrared, red, blue, and short-wave infrared bands of the image, the original time series of NDVI, EVI, and LSWI from December 1, 2022 to January 31, 2024 were calculated pixel by pixel, and denoted as the original NDVI sequence S. n EVI original sequence S e LSWI original sequence S l .

[0079] Subsequently, the missing values ​​were filled using the Spline interpolation method, and the SG filtering algorithm was introduced to smooth the interpolated sequence. The window size was 9, the polynomial order was 3, and a time series curve with a time resolution of 1 day was constructed.

[0080] Finally, by removing the data from December 2022 and January 2024 to avoid the influence of boundary effects, we obtain the 2023 remote sensing index time series dataset D={D n D e D l}, where D n To study the annual NDVI time series, D e To study the annual EVI time series, D l To study the annual LSWI time series.

[0081] Taking a rice planting pixel at latitude and longitude of 120.2776°E, 33.7332°N as an example, the time series curve of the remote sensing index for this pixel is as follows: Figure 2 As shown in the figure. Its 2023 remote sensing index time series curve is as follows. Figure 2 As shown in (a), Figure 2 In (a), the three curves represent D respectively. e Time series curve, D n Time series curves and D l Timing curve.

[0082] Step S2: Determine the effective planting season and planting system.

[0083] First, let's look at the EVI time series D for the example cells mentioned above in 2023. e The 365 vegetation indices were sorted from smallest to largest, and the 80th percentile was calculated to be 0.6799. This value was determined to be the peak threshold Th for the planting season. p .

[0084] Identifying EVI timing curves and Th p The intersection points yield 4 intersection points, with corresponding day sequences DOY 112, 143, 216, and 258. Each pair of intersection points from front to back constitutes an initial planting season, namely Initial Planting Season 1 (DOY 112~143) and Initial Planting Season 2 (DOY 216~258).

[0085] The peak EVI values ​​for the two initial planting seasons were extracted: the peak A for initial planting season 1 corresponds to day number DOY=130, with an EVI of 0.7286; the peak B for initial planting season 2 corresponds to day number DOY=233, with an EVI of 0.9576.

[0086] Verify the peak constraint conditions: ① Time constraint: The day sequence corresponding to the two peaks A and B is within the crop growth period Tg of "five consecutive days with an average daily temperature greater than 5℃", and the crop growth period Tg is DOY60~304; ② Amplitude constraint: The EVI values ​​of peaks A and B are both greater than 0.4; ③ Interval constraint: The time interval between peak A and peak B is 103 days, which is greater than 60 days.

[0087] Therefore, both initial planting seasons satisfy all constraints and are determined to be effective planting seasons. The number of effective planting seasons, m=2. Thus, the planting system for this pixel is double-cropping. Figure 2 As shown in (b).

[0088] Step S3: Extract the effective planting season start date (SOS).

[0089] First, extract the peak EVI information for two effective planting seasons: For the first effective planting season, the EVI... max-v 1 =0.7286, DOY max-v 1 =130; Second effective planting season, EVI max-v 2 =0.9576, DOY max-v 2 =233.

[0090] Secondly, extract SOS1 from the first effective planting season: using DOY max-v 1 Starting from 130, traverse the EVI timing sequence D backwards. e Find the minimum value EVI min-s 1 =0.09595; according to the formula Th s1 =0.09595+0.2×(0.7286-0.09595)=0.2225 to calculate the SOS extraction threshold; starting from DOY=130, the first day sequence with an EVI value of 0.2225 is 64, therefore SOS1=64.

[0091] Extracting SOS2 for the second effective planting season: Finding the peak day sequence DOY for the first effective planting season. max-v 1 =130, iterate through the EVI timing D within the interval DOY130~233. e Find the minimum value EVI min-s 2 =0.2035; according to the formula Th s2=0.2035+0.2×(0.9576-0.2035)=0.3543 to calculate the SOS extraction threshold; starting from DOY=233, the first day sequence with an EVI value of 0.3543 is 205, therefore SOS2=205.

[0092] Therefore, in this embodiment, the SOS set for the effective planting season is SOS. -set ={64,205}, such as Figure 2 As shown in (c).

[0093] Step S4: Identify the rice season.

[0094] First, define the time window: for SOS1=64, the time window W1=[64-20,64+20]=[44,84]; for SOS2=205, the time window W2=[205-20,205+20]=[185,225].

[0095] Then, the number of consecutive flooded days is counted: within W1, the daily determination function F(DOY) is calculated day by day. There are no cases where F(DOY)=1 is satisfied for 8 consecutive days or more, so D1=0; within W2, the number of consecutive days with F(DOY)=1 is 21, so D2=21.

[0096] Finally, the rice season is determined as follows: D1 = 0 < 8, the first valid planting season is determined to be a non-rice season; D2 = 21 ≥ 8, the second valid planting season is determined to be a rice season. The set of rice seasons is R. -set ={the second effective planting season}, n=1, such as Figure 2 As shown in (d).

[0097] Step S5: Extract key parameters for the rice growing season.

[0098] First, extract the peak EVI information for the rice season: EVI of the rice season (the second effective planting season). max-r 1 =0.9576, DOY max-r 1 =233.

[0099] Secondly, extract EOS: This rice season is the last valid planting season of the year, i.e., j=2=m, using DOY max-r 1 Starting from 233, traverse the EVI timing sequence D backwards. e Find the minimum value EVI min-rear =0.1156; according to the formula Th end =0.1156+0.2×(0.9576-0.1156)=0.2695 to calculate the EOS tail extraction threshold; starting from DOY=233, the first day sequence with an EVI value of 0.2695 is 309, therefore EOS1=309.

[0100] Then, calculate the growing season duration: T1 = EOS1 - SOS2 = 309 - 205 = 104 days;

[0101] Finally, the key parameter set for the rice growing season, Para={Para1}, was derived, where Para1={205,309,104}, as follows: Figure 2 As shown in (e).

[0102] Accuracy Verification and Result Analysis: To determine the accuracy of this method in extracting key phenological stages of rice, it was compared with traditional remote sensing methods. Traditional remote sensing methods often use thresholding based on vegetation index (EVI) time-series curves to extract key phenological stages of rice. Unlike the method of this invention, these methods cannot determine the number of effective planting seasons and do not consider the impact of complex planting systems on rice seasons and key phenological stage identification. Furthermore, unlike the method of this invention which uses bidirectional temporal detection to extract key phenological stages, this method directly obtains a threshold from the maximum and minimum values ​​of the EVI time-series curve, then judges the EVI time-series curve from front to back, identifying the first point intersecting with the threshold as SOS. This makes it highly susceptible to "interference peaks" and the influence of other crops, and the extracted SOS is easily located outside the rice growing season, leading to misjudgment of the rice season and reduced accuracy in key phenological stage identification. Again, taking pixels containing rice planting data as an example, such as... Figure 2 As shown in (f), the maximum value C of the EVI vegetation index time series curve is 0.9576, the minimum value D is 0.09595, the threshold is 0.2683, and the intersection point SOS of the threshold and the time series curve is 67. Then, taking this SOS as the starting point, it is analyzed whether there are 8 consecutive days of LSWI greater than the EVI or NDVI of the same period within the 40-day time window (47d~87d). As shown in the figure above, since it does not meet the above requirements, this pixel (non-rice + rice type) is misclassified as a single-season non-rice by traditional remote sensing methods.

[0103] The results show that the cropping systems in the farmland of the study area can be divided into two categories: single-season planting and double-season planting. Therefore, pixels containing rice cultivation can be further subdivided into two types: single-season rice cultivation and double-season crop cultivation (non-rice + rice). Pixels not containing rice cultivation can also be divided into two types: single-season non-rice cultivation and double-season crop cultivation (non-rice + non-rice). The time-series curves of the remote sensing indices for the four types are illustrated in the figure below. Figure 3 As shown, where Figure 3 (a) is a schematic diagram of the time series curves of remote sensing indices for non-rice and rice crops. Figure 3 (b) is a schematic diagram of the time series curves of the remote sensing indices for non-rice + non-rice crops. Figure 3 (c) is a schematic diagram of the time series curve of remote sensing index for single-season rice. Figure 3(d) is a schematic diagram of the time series curve of the remote sensing index for a single season of non-rice.

[0104] To quantify the accuracy of the method of this invention in identifying rice seasons under complex planting systems, 100 verification sample points in Jiangsu Province in 2023 were selected within the study area through field surveys and expert visual interpretation. Among them, 50 sample points included rice planting, specifically 25 single-season rice sample points and 25 non-rice + rice sample points; the remaining 50 sample points did not include rice planting, specifically 25 single-season non-rice sample points and 25 non-rice + non-rice (double-season non-rice) sample points.

[0105] Based on four known verification sample points, the identification accuracy of rice seasons under complex planting systems was analyzed. Specific results for the method of this invention are shown in Table 1: the identification accuracy for non-rice + rice, non-rice + non-rice, and single-season rice is 84%, the identification accuracy for single-season non-rice is 88%, and the overall classification accuracy is 85%. Traditional large-scale remote sensing methods, due to their inability to identify the number of effective crop planting seasons (i.e., the inability to determine the cropping system), have lower identification accuracy under complex planting systems. Specific results are shown in Table 2: the identification accuracy for non-rice + rice and non-rice + non-rice is 0%, the identification accuracy for single-season rice is 60%, the identification accuracy for single-season non-rice is 84%, and the overall classification accuracy is 36%. In summary, the method of this invention improves the overall classification accuracy by 49% compared to traditional remote sensing methods, indicating its applicability for identifying rice seasons under complex planting systems.

[0106] Table 1. Rice season identification accuracy of the method of the present invention under complex planting systems.

[0107]

[0108] Table 2. Accuracy of traditional remote sensing methods for identifying rice seasons under complex planting systems.

[0109]

[0110] Based on the identification of the rice season, the extraction accuracy of key phenological stages of rice was analyzed. Visual interpretation of GEE imagery yielded the DOY (Dormancy at the Beginning of Growth) time range for 30 random rice validation sample points before and after the start of rice growth. Considering the temporal resolution of the Sentinel-2 satellite (approximately 5 days), the accuracy of each method in extracting the SOS (Signal of Optimal Growth) of key phenological stages was judged by comparing its relationship with this time range. If the SOS was within the corresponding time range, the method was considered to have extracted the SOS correctly. The results showed that for the method of this invention, the SOS of 24 rice sample points was within the corresponding time range, resulting in an accuracy rate of 80% for extracting the SOS of key phenological stages. For the traditional remote sensing method, only 4 rice sample points had SOS within the corresponding time range, resulting in an accuracy rate of only 13.33%. The extraction accuracy of EOS (Earning Optical Spectrum) of key phenological stages of rice was judged similarly to that of SOS. The DOY time range for the 30 random rice validation sample points before and after the end of rice growth was obtained, and the accuracy of each method in extracting the EOS was judged by comparing its relationship with this time range. The results show that the EOS extraction accuracy of the method of this invention is 73.33%, while the EOS extraction accuracy of the traditional remote sensing method is 30%. In summary, the method of this invention improves the accuracy of SOS and EOS extraction in key phenological stages by 66.67% and 43.33% respectively compared with the traditional remote sensing method, indicating that it has a better extraction effect and can effectively reduce the impact of complex planting systems on key phenological stages of rice.

[0111] The above description is only a preferred embodiment of the present invention. It should be noted that for those skilled in the art, several improvements and modifications can be made without departing from the technical principles of the present invention, and these improvements and modifications should also be considered within the scope of protection of the present invention.

Claims

1. A method for extracting rice growth season by remote sensing based on time-series bidirectional reflectance, characterized in that, The method comprises the following steps: Step S1: constructing a remote sensing index time series dataset, specifically comprising: Step S1.1: extracting a set of farmland pixels S in the study area based on a high spatial resolution land use classification product; Step S1.2: After the spectral reflectance remote sensing image covering the study area is preprocessed, the original time series of normalized vegetation index NDVI, enhanced vegetation index EVI and land surface water index LSWI are calculated pixel by pixel, respectively denoted as NDVI original sequence S n , EVI original sequence S e , LSWI original sequence S l ; Step S1.3: Fill in the missing values in S n , S e , and S l using Spline interpolation method respectively to obtain the interpolated sequences S n ', S e ', and S l '. Step S1.4: Smooth the interpolated sequence S n ' using Savitzky-Golay filter algorithm, finally retain the valid data of the study year, form the time resolution of 1 day remote sensing index time series data set D = {D e ', S l ', S n ', S e ', S l}, wherein D n is the NDVI time series of the study year, D e is the EVI time series of the study year, and D l is the LSWI time series of the study year; Step S2: determining an effective planting season and a planting system, specifically comprising: Step S2.1: Sort all the vegetation index values in the EVI time series D e for the study year E from small to large, and take the 80% quantile as the peak threshold Th p for the growing season; Step S2.2: identifying the EVI time series D of the study year E e The intersection of the curve and Th p , and each two consecutive intersection points from front to back constitute an initial planting season S ini , and the initial planting season set is denoted as S ini-set ={S ini 1 ,S ini 2 ,...,S ini k}, and k is the number of initial planting seasons; Step S2.3: For each initial planting season S ini , extract its corresponding EVI peak, denoted as EVI peak P i , and verify it by the following three constraints: time constraint, amplitude constraint and interval constraint; Step S2.4: retain the initial planting season that meets all the above constraints, denoted as the valid planting season S val , the set of valid planting seasons is S val-set ={S val 1 ,S val 2 ,...,S val m}, m is the number of valid planting seasons, that is, the planting system of the pixel, m=1 is single-season planting, and m≥2 is multi-season planting; Step S3: extracting a starting period SOS of the effective planting season, specifically comprising: Step S3.1: For each valid planting season S val , extract its EVI maximum in the study year EVI time series D e , denoted as the planting season EVI peak EVI max-v , and the corresponding day of year as DOY max-v ; Step S3.2: for different planting systems, taking the day sequence corresponding to the EVI peak value of each effective planting season as the starting point, traversing the EVI time series forward, determining the EVI minimum value, and then calculating the extraction threshold to detect backward, the first day sequence reaching the extraction threshold is the starting period SOS of the planting season; Step S3.3: Note the set of SOSs for all valid planting seasons as SOS -set = {SOS1, SOS2,..., SOS m} ; Step S4: identifying a rice season: taking the SOS of each effective planting season as the time origin, defining a time window of 20 days before and after, counting the number of days in the window during which the LSWI is continuously higher than the NDVI or EVI at the same period, and if the number of continuous days is greater than or equal to 8, determining that the planting season is a rice season; Step S5: extracting key parameters of the rice growing season: for the rice season, according to the difference in planting system, taking the day sequence corresponding to the EVI peak value of the rice season as the starting point, traversing the EVI time series backward, determining the EVI minimum value; then calculating the tail extraction threshold, and then detecting backward to obtain the end of the rice growing season EOS, calculating the growing season length by the difference between EOS and SOS, and completing the extraction of the rice growing season.

2. The method according to claim 1, wherein, In the step S1, the calculation formulas of the vegetation index NDVI, the enhanced vegetation index EVI, and the land surface water index LSWI are as follows: ； ； ； wherein, and are near-infrared, red, blue and short-wave infrared spectral bands, respectively.

3. The method according to claim 1, wherein, In the step S2, the peak constraint condition is specifically: Time constraint: peak P i Corresponding day sequence is in the crop growing period of 5 consecutive days with daily average temperature greater than 5℃ Tg In; Amplitude constraint: peak P i Corresponding enhanced vegetation index EVI Value > 0.4; Interval constraint: if k≥2, the maximum EVI peak within the year is the main peak P max and the rest are the secondary peaks P suβ then P suβ and P max The corresponding day-of-year difference ΔDOY = |DOY suβ - DOY max | > 60 days.

4. The method according to claim 3, wherein, In the step S3.2, the SOS signal extraction rules for different planting systems are as follows: Step S3.2.1: if it is single-season planting, then it comprises: Step S3.2.1.1: Starting with DOY max-v , traverse the study year EVI time series D e forward to find the minimum value, denoted as single-season EVI minimum EVI min-s ; Step S3.2.1.2: Calculate the SOS extraction threshold Th s = EVI min-s + 0.2 x (EVI max-v - EVI min-s ) according to the formula Th s ; Step S3.2.1.3: Detecting the EVI time series D from the DOY max-v e The first point corresponding to Th s , i.e. the day order, is the start of the growing season SOS s of the effective growing season.​ Step S3.2.2: if it is multi-season planting, then it comprises: Step S3.2.2.1 : Extracting the SOS1 for the first effective planting season S val 1 , its SOS1 is extracted by the method of step S3.2.1; Step S3.2.2.2: For the t-th valid planting season S in multi-season planting val t t≥2, comprising: Step S3.2.2.2.1 : Find the t-1th valid planting season S val t-1 the EVI peak EVI max-v t-1 and the corresponding day of year DOY max-v t-1 ; Step S3.2.2.2.2: Traverse the EVI time series D of the study year e to find the minimum value of D max-v t-1 in the interval [DOY max-v t , DOY e ], denoted as the EVI minimum value EVI min-m t of the t-th valid planting season. Step S3.2.2.2.3: Calculate the t-th planting season SOS extraction threshold Th st = EVI min-m t + 0.2(EVI max-v t - EVI min-m t ), where EVI st is the t-th planting season EVI Step S3.2.2.2.4: Detecting the study year EVI time series D forward from DOY max-v t e The first point equal to Th st corresponding to the day order, i.e. the start of the growing season SOS val t of S t .​ 5. The method according to claim 4, wherein, The step S4 specifically comprises: Step S4.1: defining the start of the growing season SOS for the i-th valid planting season i , i = 1, 2,..., m; defining the time window W i from "DOY = SOS i - 20" to "DOY = SOS i + 20" Step S4.2: determining the window W of the i-th effective planting season of rice i , and comparing D l with D n , D l , and D e , respectively, and counting the number of consecutive days that satisfy D l >D n or D l >D e , and recording the number of days as the duration of continuous flooding D i of the i-th window. Step S4.3: If D i If the value is ≥8, then SOS is determined. i Corresponding effective planting season S val i For rice season R i Conversely, it is a non-rice season; let the set of rice seasons be R. -set ={R1,R2,...,R n }, where n is the number of rice seasons, and n≤m.

6. The method according to claim 5, wherein, The step S5 specifically comprises: Step S5.1: Extract the EVI peak value of the jth rice season R j , j = 1, 2, …, n, in the EVI time series D e of the study year E, denoted as the rice season EVI peak value EVI max-r j , and the corresponding day sequence is denoted as DOY max-r j ; Step S5.2: If R j is the last valid rice planting season of the year, i.e. j = m, then includes: Step S5.2.1: Take DOY max-r j as the starting point, traverse back to find the minimum value of the study year EVI time series D e , recorded as the minimum value of the end of the rice season EVI, EVI min-rear ; Step S5.2.2: Calculate the EOS tail extraction threshold Th end = EVI min-rear + 0.2 x (EVI max-r j - EVI min-rear ) according to the formula Th end ; Step S5.2.3: Detecting EVI time series D backward from DOY max- r j as the starting point e , the first point corresponding to the day sequence equal to Th end , that is, the growth season end period EOS j of R j 1 ; Step S5.3: If R j If there is a valid planting season later, i.e. j < m, then include: Step S5.3.1 : Find R j the next valid growing season S val j+1 the EVI peak EVI max-v j+1 and the corresponding day of year DOY max-v j+1 ; Step S5.3.2: traverse the EVI time series D of the study year E e to find the minimum value of D max-r j in the interval [DOY max-v j+1 , DOY e ], denoted as the middle EVI minimum value EVI min-rmid of the rice season. Step S5.3.3: According to formula Th end j =EVI min-rmid +0.2×(EVI max-r j -EVI min-rmid ), calculate the EOS tail extraction threshold Th for the j-th rice season. end j ; Step S5.3.4: Detect D max-r j with DOY e as starting point backwards end j , the first point equal to Th j corresponds to the end of season EOS j of R Step S5.4: Uniform record R j EOS j , j = 1, 2,..., n, R j is calculated according to the formula j = EOS j - SOS j ; Step S5.5: integrate all R j of SOS j , EOS j , T j to form a set of key parameters of rice growing season Para={Para1, Para2,..., Para n}, wherein Para j ={SOS j , EOS j , T j}, and the extraction of rice growing season is completed.

7. The method according to claim 6, wherein, In the step S3.2, the forward detection research year EVI time series D e The specific way is: if it is single season planting, from DOY max-v , read the value of D e in the order of day sequence from day to day, when the read value is equal to Th s for the first time, stop detection, and record the day sequence corresponding to the value. If it is multi-season planting, from DOY max-v t , read the value of D e every day in the direction of decreasing day order, and when the read value is first equal to Th st , stop detection and record the day order corresponding to the value.

8. The method according to claim 7, wherein, The step S4.2, the continuous flooding days D i The statistical rule is: within W i , when the condition of D l >D n or D l >D e occurs, the time is started to be counted, if the subsequent date continuously satisfies the condition, the continuous counting is continued, if a date does not satisfy, the counting is reset, and finally the longest continuous counting days within W i are taken as D i .

9. The method according to claim 8, wherein, The step S5.2.3 and step S5.3.4, backward detection of the EVI time series D of the study year E e The specific way is: from DOY max-r j , read the value of D e day by day in the direction of increasing day order, when the read value is equal to Th end or Th end j for the first time, stop detection, and record the day order corresponding to the value.

Images