Mapping Ratoon Rice Fields Based on SAR Time Series and Phenology Data in Cloudy Regions

Abstract

Ratoon rice (RR) has emerged as an active adaptation to climate uncertainty, stabilizing total paddy rice yield and effectively reducing agriculture-related ecological environmental issues. However, identifying key remote sensing parameters for RR under cloudy and foggy conditions is challenging, and existing RR monitoring methods in these regions face significant uncertainties. Here, given the sensitivity of synthetic aperture radar (SAR) backscattering signals to the crop phenological period, this paper introduces a threshold model utilizing Sentinel-1A SAR data and phenological information for mapping RR. The Yongchuan District of Chongqing, which is often cloudy and foggy, was selected as a specific study region where VH-polarized backscatter coefficients of Sentinel-1 images were obtained at 10 m spatial resolution in 2020. Based on the proposed threshold model, the RR extraction overall accuracy was up to 90.24%, F1 score was 0.92, and Kappa coefficient was 0.80. Further analysis showed that the extracted RR boundaries exhibited high consistency with true Sentinel-2 remote sensing images and the RR extracted area was in good agreement with the actual planted area situation. This threshold model demonstrated good applicability in the studied cloudy and foggy region, and successfully distinguished RR from other paddy rice types. The methodological framework established in this study provides a basis for extensive application in China and other significant RR-producing regions globally.

1. Introduction

Ratoon rice (RR) is a mode of paddy rice cultivation. After the first rice crop is harvested, the surviving dormant buds on the stubble are used for reproduction, ultimately achieving two harvests in a given period. Through specific cultivation management measures, the buds are induced to sprout again, tiller, panicle, mature, and be harvested, thus completing the paddy rice growth cycle. RR possesses the advantages of shorter growth duration, higher yield, and water and labor savings compared with traditional early rice, late rice, and middle rice cropping systems [1,2,3,4]. In recent years, RR has been vigorously promoted in south-western China, with the cultivation area expanding steadily since 2011 [5]. It has emerged as a new planting system in the fight to adapt to climate change, increase planting frequency, stabilize total paddy rice production, and mitigate agricultural environmental issues [6].

Field investigation is the traditional method by which the area and planting regions of paddy rice are obtained [7]. However, this method is not feasible for large-scale studies because it is time consuming and labor intensive. The challenge in promptly acquiring accurate data on the geographical distribution of paddy rice fields stems from the labor-intensive nature of conducting field surveys. The emergence of remote sensing technology has addressed this issue by offering advantages such as rapid information acquisition, wide detection coverage, and relatively low cost. Remote sensing applications in the agricultural sector have become increasingly common, providing the most effective method for extracting the planting area of RR. This technology holds great importance for the large-scale continuous management of agricultural production.

Optical remote sensing (ORS) and synthetic aperture radar (SAR) data are commonly employed in the mapping of both double-season and single-season rice, yet their utilization in RR mapping is comparatively limited. Different vegetation parameters can be obtained from ORS data [8], including the normalized difference vegetation index, enhanced vegetation index, and modified normalized difference water index, among others. These parameters, derived from MODIS [9], Landsat [10], and Sentinel-2 [11] data, have been effectively employed in the delineation of paddy rice fields. Currently, ORS is the mainstream data source for monitoring RR. For example, Liu et al. [12] and Li et al. [13] combined phenological information with Sentinel-2 or MOD09Q1 images to map RR in Hubei Province, achieving a high overall accuracy (OA) (Sentinel-2: 0.76, MOD09Q1: 0.94) and Kappa coefficient (Sentinel-2: 0.65, MOD09Q1: 0.81). Chen et al. [14] introduced a novel index by integrating Landsat and Sentinel-2 data to map RR in Qichun County, Hubei Province, achieving an OA of 0.80. Similarly, in our previous study, a fused time-series optical dataset and a threshold model were applied to monitor the current status of RR growth in Yongchuan District, Chongqing, and the results showed that the overall accuracy (90.73%) and Kappa coefficient (0.81) were at a high level [15]. However, a substantial drawback of ORS data lies in its susceptibility to interference from clouds and fog. Consequently, the availability of such data may be limited, particularly in regions characterized by persistent cloud cover, such as southwestern China [16]. In contrast, SAR possesses the capability of 24-h operation, regardless of weather conditions, penetrating cloud cover effectively, and exhibiting minimal sensitivity to meteorological fluctuations [17]. There has been an increasing trend in utilizing SAR data for monitoring paddy rice fields over recent years because SAR backscatter time series can be used as a basis for identifying key features of the paddy rice growth period. As early as 1996, SAR data from the European Remote Sensing Satellite 1 (ERS-1) was utilized in the monitoring of paddy rice [18]. However, early research was limited by the scarcity of high-quality real images, limiting the scope to localized regions. With the advent of Sentinel-1 SAR data, offering a spatial resolution of 10 m and revisiting an area every 12 days per satellite, the utilization of this resource now holds great promise for precisely mapping various types of paddy rice.

The utilization of Sentinel-1 SAR data has been employed to map paddy rice in diverse regions and cropping systems, mainly encompassing double-season and single-season rice cultivation [8,19,20,21]. The prevalent approach in paddy rice mapping involves developing a threshold model or employing machine-learning techniques. The principle of a threshold-based decision tree is to establish mathematical formulae based on the dynamic range of radar backscattering or variance during the paddy rice growth period. The threshold model has a wide range of applications and can achieve >80% accuracy [8,15,22]. Computer technology advancements have led to the incorporation of machine-learning techniques in SAR data monitoring, facilitating the extraction of distinct values from changes in backscatter coefficients throughout the growth stages of paddy rice. These values serve as inputs for different machine-learning techniques, including the decision tree [23,24], support vector machine [25], and random forest [26,27,28] models; these models yield paddy rice monitoring accuracies exceeding 80%. Hence, both threshold models and machine-learning methods can achieve effective monitoring of different regions for single-season and double-season rice cultivation. The expense of training machine-learning models can be significant, requiring substantial amounts of training data to yield desirable outcomes [29] Moreover, there remains a scarcity of practical research concerning remote sensing monitoring of RR. In contrast, the idea of threshold models is relatively simple; once the thresholds derived from the time-series SAR data have been set, the process of identifying paddy rice within the study region can commence [29].

To enhance the monitoring system for RR across diverse paddy rice types and overcome the challenges of ORS in cloudy and foggy regions, based on our previous research [15], we set out to achieve the following goals in this study: (1) We selected the VV and VH polarizations of Sentinel-1A SAR data, then applied the Lee sigma filter to control speckle noise in the SAR data. Subsequently, we employed the Savitzky–Golay (S-G) filter to temporally smooth the SAR image sequences. (2) To develop a threshold-based model for monitoring RR in regions characterized by cloudy and foggy weather conditions. By analyzing the temporal patterns of backscatter coefficients for diverse land-cover types and combining them with the phenological calendar, we determined key windows and thresholds to extract RR in the study region.

2. Study Region and Data

2.1. Study Region

Yongchuan District lies in the western part of Chongqing, China, situated on the northern bank of the upper Yangtze River. Its geographic coordinates range from approximately 28°56′ to 29°34′N and 105°38′ to 106°05′E, encompassing a total area of 1576 km² (Figure 1a). Yongchuan District is primarily characterized by landforms such as low mountains, hills, and terraces. The district experiences a subtropical monsoon humid climate, with average annual precipitation of 1015.0 mm, and an average temperature of 17.7 °C. Paddy rice is the primary crop cultivated in Yongchuan District, in a manner typical of paddy rice cultivation in mountainous and hilly regions. According to the Chongqing Statistical Yearbook for 2021 (https://tjj.cq.gov.cn/zwgk_233/tjnj/tjnj.html?url=http://tjj.cq.gov.cn/zwgk_233/tjnj/2021/indexch.htm, accessed on 15 December 2021), Yongchuan District ranks fifth in grain production in Chongqing. The average annual paddy rice planting area is 400 km², constituting approximately 60% of the total cereal crop planting area. The total paddy rice production is approximately 3.4 × 10⁸ kg, representing approximately 70% of the total grain production in Yongchuan District (http://www.cqyc.gov.cn/zwgk_204/zfxxgk/zfxxgkml/tjxx/, accessed on 20 June 2023). According to historical weather forecast data from the China Meteorological Network (https://www.cma.gov.cn/, accessed on 9 June 2021), the weather conditions experienced in Yongchuan District during 2020 were as follows: overcast weather accounted for 41%, partly cloudy weather accounted for 35%, rainy weather accounted for 9%, and sunny days accounted for only 15% (Figure 1c). These data indicate that Yongchuan District experienced cloudy or rainy weather for 85% of the year, classifying it as a typical region with frequent cloud cover and rainfall. Therefore, selecting this region for RR remote sensing monitoring research holds significance in promoting and demonstrating RR mapping capabilities in cloudy and foggy regions.

2.2. Data

The data utilized here mainly included the 2020 Sentinel-1A time-series data of Yongchuan District, field survey data, statistical data on paddy rice planting area, paddy rice planting calendar, and relevant auxiliary land-cover data.

2.2.1. Sentinel-1 Data

The Sentinel-1 mission comprises a pair of satellites, Sentinel-1A and Sentinel-1B, and plays a vital role in the European Space Agency’s Copernicus Programme for Earth observation. Sentinel-1 is outfitted with C-band SAR, allowing for continuous imaging under all weather conditions at any time of the day. The revisit cycle for each satellite is set at 12 days, and theoretically, the revisit cycle for both satellites can be as short as six days. Sentinel-1 SAR offers four operational modes: interferometric wide swath (IW), extra wide swath (EW), wave (WV), and stripmap (SM). The polarization modes of Sentinel-1 include single and dual polarization. Single polarization options are either HH or VV, while dual polarization options are represented by either HH + HV or VV + VH. We accessed Sentinel-1A SAR data through the Google Earth Engine (GEE) platform, at 10 m spatial resolution; these data had already been preprocessed, including orbit correction, thermal noise removal, radiometric calibration, and terrain correction. For this study, the acquisition mode was IW mode, and the product type was ground range detected imagery, including both VV and VH polarizations. Sentinel-1B data were missing over the period February–November 2020 in the study region. Therefore, 24 Sentinel-1A images were ultimately acquired, with a revisit period of 12 days.

2.2.2. Field Data Collection

Field survey data were used to construct the threshold model to monitor RR and to validate the accuracy of the RR mapping result. We utilized a real-time kinematic surveying device for our first field survey on 31 May and 1 June 2021, collecting samples of both rice and non-rice points. Our second field survey was conducted on 29 and 30 August 2021, to collect sample points of single-season rice and RR. We selected large and flat fields for sampling to match the scale of remote sensing data. The center of each field was used as the position of each sampling point, which was projected to Sentinel-1 SAR data. Figure 1a shows the 185 RR sample points and 76 non-RR sample points (including water body, forest, urban, and others). The training and testing sets of RR samples were randomly divided into a ratio of 3:7. Non-RR sample points do not need to be divided into training and test sets. In the construction of threshold models, they are used to mask background land-cover information. Through our field investigations, we identified some issues with RR cultivation in Yongchuan District (Figure 1b), including discontinuous planting of RR and damage to some rice stubbles during the harvest stage of the main rice crop, posing challenges to the remote sensing monitoring of RR. However, overall, RR cultivation in the south of Yongchuan District was better than that in the north, with mostly contiguous planting regions in the former. The majority of the RR exhibited robust growth, characterized by full stubble coverage, dense regrowth buds, and a notable abundance of paddy rice panicles (Figure 1b), indicating that Yongchuan District benefits from mature RR cultivation management techniques.

2.2.3. RR Planting Calendar and Other Ancillary Data

The RR planting calendar for Yongchuan District, as shown in Figure 2, was obtained from the Chongqing General Agricultural Technology Station (http://www.cqates.com/, accessed on 20 March 2021). In this region, farmers begin seedling cultivation in March, followed by land leveling, field irrigation to a depth of 2–15 cm, and fertilization, preparing for transplanting one month later. In April, during the flooding stage of the first rice crop, the water depth is typically maintained at approximately 3 cm. As the first rice crop enters the tillering stage (May), followed by the heading and filling stages (June–July), the biomass of the crop increases, requiring more water, with the water depth increasing to 5–10 cm [19]. Timely drainage before maturity is necessary to prepare for harvesting (mid-August). During the harvest of the first rice crop, it is recommended to leave stubble of 20–40 cm [30]. Subsequent fertilization and shallow water irrigation in late August to early September can help stimulate the dormant buds to sprout and promote growth. Over approximately 60 days, the ratoon crop goes through the stages of tillering, heading, filling, maturity, and harvesting. RR grows rapidly, with its growth period being considerably shorter than that of the first rice crop.

Other ancillary data included the planting area of RR, land-cover data in 2020, and Google Earth high-resolution images of Yongchuan District. The planting area of RR was obtained from official government documents of Yongchuan District and professional researchers, and it was used to evaluate the accuracy of RR monitoring results. We utilized the high-spatial-resolution (10 m) land-cover product of ESRI generated from Sentinel-2 images, and the OA of this dataset was 0.85 [31]. This product was employed to extract land-cover masks for urban land, water body, forest, and grassland in the study region. High-resolution images from Google Earth were used to assist in verifying the accuracy of the ground sample points.

3. Methodology

The main schematic flow of RR mapping is illustrated in Figure 3. First, we performed Lee sigma filtering on Sentinel-1A SAR data to mitigate speckle noise. Then, we applied S-G filtering to smooth the Sentinel-1A time-series data of the study region. This paper compared the time-series curves of RR using the VV and VH polarization modes of Sentinel-1A SAR data and selected the optimal polarization mode. Subsequently, we integrated the RR phenological information and 30% of the RR training sample points for Sentinel-1A time-series feature analysis. Based on analysis results, we then identified the characteristics of RR and extracted relevant thresholds to preliminarily determine the distribution of RR. Finally, we utilized land-cover data as a mask to filter out potential sources of interference, acquiring the final distribution map of RR cultivation, before conducting accuracy verification.

3.1. Preprocessing and Polarization Mode Selection

Given that RR plants across the study region will likely exhibit variations in phenological features, we calculated the standard deviation and used the error band diagram to reflect the error distribution range of SAR intensity around the average value (Figure 4). The Lee sigma filter is effective in removing speckle noise from Sentinel-1A data [8,32,33]; we used SNAP 8.0 software for Lee sigma filtering in this study, as depicted in Figure 4. However, despite this processing, the curve still exhibited saw-toothed interference patterns. This suggests that employing only Lee sigma filtering is insufficient for obtaining RR growth time-series curves with smooth and characteristic features. Hence, further smoothing using S-G filtering was conducted. Figure 4 shows that S-G filtering method effectively mitigated interference while preserving important curve characteristics, such as prominent peaks and troughs [32,34].

The backscatter time-series curves of RR exhibited similar characteristics in both VH and VV polarization modes (Figure 4): the first trough occurred during the transplanting stage, while the second trough appeared during the stubble stage. However, compared with VV polarization, VH polarization exhibited more pronounced curve characteristics, with a greater overall variation in the curve, indicating that the variations in backscatter response under VH polarization exhibit greater sensitivity to the development of RR [8,35,36]. Therefore, we employed the backscatter coefficient under VH polarization mode to construct our threshold model for RR.

3.2. RR Monitoring Threshold Model

The primary time windows of SAR data were concentrated around the transplanting, heading, filling, and harvesting stages of the first rice crop, and the stubble, heading, and filling stages of the ratoon crop. Similar to Figure 4, we calculated the standard deviation and used the error band diagram to reflect the error distribution range of SAR intensity around the average value of the five land-cover types (Figure 5). As shown in Figure 5, during the transplanting stage, the occurrence of a strong flooding event resulted in a reduction in the VH backscatter coefficient, approximating the backscatter coefficient of water body. This phenomenon manifested as the first ‘V’-shaped trough feature. As the paddy rice continued to grow, there was an increase in the VH backscatter coefficient. For the tillering and heading periods of the first rice crop, the VH backscatter coefficient approached that of vegetation, resulting in the first peak feature. In the maturity and harvesting stages of the first rice crop, the VH backscatter coefficient gradually decreased, exhibiting a secondary ‘V’-shaped trough feature. The retained paddy rice stubble contains certain biomass information; hence, the backscatter coefficient at this stage was higher than that of the first ‘V’-shaped trough. For the heading and filling stages of the ratoon crop, there was a notable increase in the VH backscatter coefficient, culminating in the emergence of the second peak feature. As the RR matured and was harvested, there was a reduction in the VH backscatter coefficient. In general, the VH backscatter curve exhibited two troughs and two peaks during the growth period of RR. Urban and forest land-cover types have higher backscatter values, while water body demonstrate lower backscatter coefficients, making them easily distinguishable from paddy rice and RR. On the basis of the key growth stages and phenological calendar of RR, we designed an SAR monitoring threshold model (Figure 6): Sentinel-1A VH time series data were constructed by combining sentinel-1A images and RR training samples. The thresholds of different levels were determined by experiments. After satisfying the threshold conditions, the RR was gradually classified, and those that did not meet the threshold conditions were classified as non-RR. The steps were as follows:

(1) Masking background land-cover information. The backscatter coefficient of water body is lower than those of urban and forest land [8], while the backscatter coefficient of RR falls within these upper and lower limits. We calculated the average backscatter coefficient of water body as the minimum value and the average backscatter coefficient of urban land as the maximum value. If the average value of a pixel was not within the range of the maximum and minimum values, then this pixel was considered to be non-RR; otherwise, we proceeded to the next step:

$$\overline{{\mathrm{V}\mathrm{H}}_{\mathrm{y}\mathrm{e}\mathrm{a}\mathrm{r}}}>\mathrm{T}1 \mathrm{or} \overline{{\mathrm{V}\mathrm{H}}_{\mathrm{y}\mathrm{e}\mathrm{a}\mathrm{r}}}<\mathrm{T}2$$

(1)

here, $\overline{{\mathrm{V}\mathrm{H}}_{\mathrm{y}\mathrm{e}\mathrm{a}\mathrm{r}}}$ represents the average VH backscatter coefficient of a pixel throughout the year; $\mathrm{T}1$ is the average backscatter coefficient of water body, and $\mathrm{T}2$ is the average backscatter coefficient of urban land.

(2) Identifying two ‘V’-shaped trough features. During the growth stage of RR, the backscatter coefficient curve exhibited two trough values. We utilized the slopes on both sides of these ‘V’-shaped trough features to identify RR:

$${\mathrm{S}\mathrm{l}\mathrm{o}\mathrm{p}\mathrm{e}}_{\mathrm{f}\mathrm{i}\mathrm{r}\mathrm{s}\mathrm{t}–\mathrm{l}\mathrm{e}\mathrm{f}\mathrm{t}}<\mathrm{T}3 \mathrm{a}\mathrm{n}\mathrm{d} {\mathrm{S}\mathrm{l}\mathrm{o}\mathrm{p}\mathrm{e}}_{\mathrm{f}\mathrm{i}\mathrm{r}\mathrm{s}\mathrm{t}–\mathrm{r}\mathrm{i}\mathrm{g}\mathrm{h}\mathrm{t}}<\mathrm{T}4$$

(2)

$${\mathrm{S}\mathrm{l}\mathrm{o}\mathrm{p}\mathrm{e}}_{\mathrm{s}\mathrm{e}\mathrm{c}\mathrm{o}\mathrm{n}\mathrm{d}–\mathrm{l}\mathrm{e}\mathrm{f}\mathrm{t}}<\mathrm{T}5 \mathrm{a}\mathrm{n}\mathrm{d} {\mathrm{S}\mathrm{l}\mathrm{o}\mathrm{p}\mathrm{e}}_{\mathrm{s}\mathrm{e}\mathrm{c}\mathrm{o}\mathrm{n}\mathrm{d}–\mathrm{r}\mathrm{i}\mathrm{g}\mathrm{h}\mathrm{t}}<\mathrm{T}6$$

(3)

where, ${\mathrm{S}\mathrm{l}\mathrm{o}\mathrm{p}\mathrm{e}}_{\mathrm{f}\mathrm{i}\mathrm{r}\mathrm{s}\mathrm{t}–\mathrm{l}\mathrm{e}\mathrm{f}\mathrm{t}}$ and ${\mathrm{S}\mathrm{l}\mathrm{o}\mathrm{p}\mathrm{e}}_{\mathrm{f}\mathrm{i}\mathrm{r}\mathrm{s}\mathrm{t}–\mathrm{r}\mathrm{i}\mathrm{g}\mathrm{h}\mathrm{t}}$ represent the slopes before and after the first ‘V’-shaped trough feature, respectively, which corresponds to the transplanting period; similarly, ${\mathrm{S}\mathrm{l}\mathrm{o}\mathrm{p}\mathrm{e}}_{\mathrm{s}\mathrm{e}\mathrm{c}\mathrm{o}\mathrm{n}\mathrm{d}–\mathrm{l}\mathrm{e}\mathrm{f}\mathrm{t}}$ and ${\mathrm{S}\mathrm{l}\mathrm{o}\mathrm{p}\mathrm{e}}_{\mathrm{s}\mathrm{e}\mathrm{c}\mathrm{o}\mathrm{n}\mathrm{d}–\mathrm{r}\mathrm{i}\mathrm{g}\mathrm{h}\mathrm{t}}$ represent the slopes before and after the second ‘V’-shaped trough feature, respectively, which corresponds to the maturity and harvesting stages of the first rice crop; and T3, T4, T5, and T6 represent undetermined threshold slopes. We traversed two time intervals on both sides of the trough values, using day of year (DOY) as the independent variable and the VH backscatter coefficient as the dependent variable, to calculate the slopes. If the curve exhibited two ‘V’-shaped features and the slope values fell within the specified threshold, then we proceeded to the next step; otherwise, the pixel was classified as non-RR.

(3) Peak and trough interval days. Peaks and troughs were extracted throughout the entire time range, and their corresponding DOY values were output for subsequent interval calculations.

$$\mathrm{T}7\le \mathrm{L}\mathrm{B}\mathrm{T}\le \mathrm{T}8$$

(4)

$$\mathrm{T}9\le \mathrm{L}\mathrm{B}\mathrm{T}\le \mathrm{T}10$$

(5)

where, T7, T8, T9, and T10 were determined based on 30% of the training samples; T7 and T8 represent the shortest and longest interval days between the two troughs, respectively, while T9 and T10 represent the shortest and longest interval days between the two peaks, respectively. These values were determined based on the mean and variance of the samples:

$$\mathrm{T}7=\overline{{\mathrm{L}\mathrm{T}\mathrm{P}}_{\mathrm{t}\mathrm{e}\mathrm{s}\mathrm{t}}}-\mathrm{i}{\mathsf{\sigma }}_{{\mathrm{L}\mathrm{T}\mathrm{P}}_{\mathrm{t}\mathrm{e}\mathrm{s}\mathrm{t}}}$$

(6)

$$\mathrm{T}8=\overline{{\mathrm{L}\mathrm{T}\mathrm{P}}_{\mathrm{t}\mathrm{e}\mathrm{s}\mathrm{t}}}+\mathrm{i}{\mathsf{\sigma }}_{{\mathrm{L}\mathrm{T}\mathrm{P}}_{\mathrm{t}\mathrm{e}\mathrm{s}\mathrm{t}}}$$

(7)

$$\mathrm{T}9=\overline{{\mathrm{L}\mathrm{B}\mathrm{P}}_{\mathrm{t}\mathrm{e}\mathrm{s}\mathrm{t}}}-\mathrm{i}{\mathsf{\sigma }}_{{\mathrm{L}\mathrm{B}\mathrm{P}}_{\mathrm{t}\mathrm{e}\mathrm{s}\mathrm{t}}}$$

(8)

$$\mathrm{T}10=\overline{{\mathrm{L}\mathrm{B}\mathrm{P}}_{\mathrm{t}\mathrm{e}\mathrm{s}\mathrm{t}}}+\mathrm{i}{\mathsf{\sigma }}_{{\mathrm{L}\mathrm{B}\mathrm{P}}_{\mathrm{t}\mathrm{e}\mathrm{s}\mathrm{t}}}$$

(9)

where, $\overline{{\mathrm{L}\mathrm{B}\mathrm{P}}_{\mathrm{t}\mathrm{e}\mathrm{s}\mathrm{t}}}$ is the average interval days between two peaks of the training samples, $\overline{{\mathrm{L}\mathrm{T}\mathrm{P}}_{\mathrm{t}\mathrm{e}\mathrm{s}\mathrm{t}}}$ represents the average interval days between two troughs of the training samples, $\mathrm{i}{\mathsf{\sigma }}_{{\mathrm{L}\mathrm{B}\mathrm{P}}_{\mathrm{t}\mathrm{e}\mathrm{s}\mathrm{t}}}$ denotes the variance of the interval days between two peaks of the training samples, and $\mathrm{i}{\mathsf{\sigma }}_{{\mathrm{L}\mathrm{T}\mathrm{P}}_{\mathrm{t}\mathrm{e}\mathrm{s}\mathrm{t}}}$ denotes the variance of the interval days between two troughs of the training samples. Here, $i$ = 1, 1.25, …, 3 times the standard deviation. Ultimately, T7, T8, T9, and T10 each had nine types of threshold result. Through mathematical combinatorial methods, 6561 combinations for T7, T8, T9, and T10 (9⁴ = 6561) were derived, culminating in the identification of the most effective threshold combination. To minimize variability in the outcomes, the training and test samples were randomly divided 100 times to obtain a more scientific threshold. T1, T2, T3, T4, T5, and T6 were also calculated using 30% of the training samples. The application process of threshold model was realized by MATLAB R2020b software. The optimal thresholds were as follows:

T1 = −9.79, T2 = −27.1593, T3 = 0.01, T4 = 0.009, T5 = 0.006, T6 = 0.01

(10)

T7 = 90 days, T8 = 156 days, T9 = 60 days, T10 = 110 days

(11)

3.3. Accuracy Metrics

The confusion matrix for RR mapping was computed using 70% of the RR testing samples and all non-RR samples. In Yongchuan District, there were a total of 205 samples for evaluation, comprising 129 RR samples and 76 non-RR samples. The accuracy of RR mapping was evaluated using OA, user’s accuracy (UA), producer’s accuracy (PA), F1 score, and the Kappa coefficient. The F1 score comprehensively reflected the harmonic mean of UA and PA, providing a thorough assessment of the result [37]. These accuracy metrics are represented by the following equations:

$$\mathrm{O}\mathrm{A}={\displaystyle \frac{S_d}{N}}$$

(12)

$$\mathrm{U}\mathrm{A}={\displaystyle \frac{S_{ij}}{S_{*j}}}$$

(13)

$$\mathrm{P}\mathrm{A}={\displaystyle \frac{S_{ij}}{S_{i*}}}$$

(14)

$$\mathrm{F}1 \mathrm{s}\mathrm{c}\mathrm{o}\mathrm{r}\mathrm{e}={\displaystyle \frac{2\times UA\times PA}{UA+PA}}$$

(15)

$$\mathrm{K}\mathrm{a}\mathrm{p}\mathrm{p}\mathrm{a}={\displaystyle \frac{OA-P_e}{1-P_e}}, P_e={\displaystyle \frac{\sum (S_{*j}\times S_{i*})}{N^2}}$$

(16)

where, $S_d$ denotes the number of samples correctly classified, $N$ represents the number of testing samples, $S_{ij}$ denotes the number of samples classified as land-cover type$j$ when the true land-cover type is $i$, $S_{*j}$ is the total number of samples with true land type $j$, and $S_{i*}$ represents the total number of samples classified as land-cover type $i$.

4. Results

4.1. RR Mapping Results and Accuracy Verification

An RR map was produced for the study region utilizing VH time-series data and the threshold model (Figure 7A). Overall, RR cultivation was extensive in Yongchuan District, primarily distributed along riverbanks and regions of flat terrain. A spatial resolution of 10 m can effectively identify fragmented fields in the southwest region, reduce the impact of mixed pixels, and enhance the reliability of the results. The actual planting area of RR in Yongchuan District ranges between 158 and 180 km² [15]. The total monitoring area for RR fields was 199.79 km², which closely aligns with the available statistical data. Figure 7B illustrates that Zhutuo Town had the largest area of RR (17.27 km²), followed by Laisu Town (17.37 km²), Sanjiao Town (14.42 km²) and Hegeng Town (14.15 km²). In addition, the townships with RR cultivation areas exceeding 10 km² included Xianlong Town, Chenshi Street, and Da’an Street. Cultivation of RR in these regions was extensive, with the total RR area in these townships accounting for 51.01% of the total RR cultivation area in Yongchuan District.

We conducted an accuracy validation for the RR mapping result using the RR testing samples and all non-RR samples (

Table 1. Performance evaluation of RR mapping based on a confusion matrix: overall accuracy (OA), producer’s accuracy (PA), user’s accuracy (UA), F1 score, and Kappa coefficient.

Item	RR (Actual Value)	Non-RR (Actual Value)	Total
RR (predicted value)	110	0	110
Non-RR (predicted value)	19	76	95
Total	129	76
OA (%)	90.24
PA (%)	85
UA (%)	100
F1 score	0.92
Kappa coefficient	0.8

), which gave an OA of 90.24% and Kappa coefficient of 0.8 (0.61 $\le$ Kappa $\le$ 0.8 indicates a high level of heterogeneity). Moreover, the PA and UA were 85% and 100%, respectively. UA represents the proportion of RR (or non-RR) samples that are correctly classified and predicted as RR (or non-RR) samples. PA represents the proportion of RR (or non-RR) samples that are correctly classified relative to the true RR (or non-RR) samples. Table 1 also shows the proposed method can get a promising result in the study region (e.g., F1 score = 0.92). These results indicate that in cloudy and rainy regions, the SAR signal and setting of threshold model parameters can effectively differentiate between RR and other land-cover types.

4.2. Detailed Spatial Features of RR

To analyze the spatial details of the RR extraction results, two random sub-areas were selected for comparison with real remote sensing images (Figure 7(a-1–a-4)). Generally, RR distributions showed good agreement with Sentinel-2 optical images. This indicates that the method outlined herein successfully differentiated between RR and non-RR land cover. However, the spatial distributions of the extracted patches were somewhat fragmented, which can be attributed to three possible factors. First, the pixel-based method for identifying RR inevitably encounters a “salt and pepper” effect. Second, radar data is susceptible to external interference, leading to some data instability; specifically, this may be related to the imaging principles of SAR, where factors such as electromagnetic environmental interference and the characteristics of ground targets can cause variations in backscatter coefficient characteristics [8,38,39]. Third, the cultivation scenario of RR in the research region was rather complex; for example, through field investigations, we found that there were differences in the growth conditions of RR in a single field plot (Figure 1b). For future improvement, an object-oriented extraction method could be employed. This method involves image segmentation clustering before extraction, considering neighborhood pixel features, which effectively enhances classification accuracy [40].

4.3. Comparison with Existing Methods

We compared the proposed threshold model with other methods of RR mapping, considering the type of method, phenological information required, data sources, spatial resolution, OA, UA, PA, F1 score, and Kappa coefficient (

Table 2. Comparison of the accuracy of four methods of RR mapping: type of method, phenological information required, data sources, overall accuracy (OA), user’s accuracy (UA), producer’s accuracy (PA), F1 score, and Kappa coefficient.

	This Study	Zhao et al. [15]	Liu et al. [12]	Li et al. [13]
Type of method	Threshold model	Threshold model	Index method	RR index
Phenological information	Yes	Yes	Yes	Yes
Data sources	Sentinel-1A; ESRI land cover	Sentinel-2; Landsat-8 OLI; MOD09GA	Sentinel-2; FROM-GLC10; DEM	MOD09Q1; MCD12Q1; DEM
Spatial resolution (m)	10	10	10	250
OA (%)	90.24	90.73	No comparison	No comparison
UA (%)	100	/
PA (%)	85	/
F1 score	0.92	/
Kappa coefficient	0.80	0.81

). This study and that of Zhao et al. [15] both used rules and thresholds to map RR. Meanwhile, during the RR growing period, index-based methods were used by Liu et al. [12] and Li et al. [13] for RR identification. The index method is a technique used to quantify and evaluate specific phenomena, characteristics, or conditions by calculating and utilizing various indices. These indices are typically quantitative indicators derived from raw data and are standardized or normalized through mathematical formulas to facilitate comparisons across different datasets, time periods, or spatial regions. In the remote sensing monitoring of RR, relevant information is extracted from raw data, and a new index is constructed through recombination to extract RR. On the other hand, a threshold, also known as a critical value, refers to the minimum or maximum value at which an effect can occur. In mathematical or statistical modeling, a threshold model is any model where a threshold value, or set of threshold values, is used to distinguish ranges of values where the behavior predicted by the model varies in some important way. These studies all require RR phenological information. Herein, we applied the dynamic range of Sentinel-1A SAR backscatter as the key feature for RR mapping. However, the other three aforementioned methods used optical images (Sentinel-2, Landsat-8 OLI, MOD09GA, and MOD09Q1), which are susceptible to clouds and fog, as their main data source. Since these models were developed and tested using different ground-truth datasets obtained in different locations, this study did not compare the accuracy evaluation indicators with those of Liu et al. [12] and Li et al. [13] (Table 2). However, at equivalent spatial resolution (10 m), same conditions of location (Yongchuan District, Chongqing), time (2020), and sample points, we carried out a comparison with the evaluation indicators of our previous study [15]. Table 2 shows that the OA and Kappa coefficients of this study are slightly lower than Zhao et al. [15], but the gap is within 0.01. The results indicate that the use of SAR data to identify RR can effectively overcome the influence of cloudy weather, but the recognition accuracy may be slightly inferior to ORS data.

5. Discussion

5.1. Performance of the Proposed Model in Classifying Different Paddy Rice Types

The main focus of this model is to identify “double peak” and “double trough” features to determine whether the analyzed component is RR. The time-series curve of single-season rice exhibits only one peak and one trough, whereas that of RR has “double peaks” and “double troughs”, allowing the two to be clearly distinguished. Double-season rice time-series curves also exhibit “double peaks” and “double troughs”, but the second-season rice of double-season rice requires re-transplantation, accompanied by a strong flooding signal. In contrast, the ratoon crop only needs to retain 20–40 cm of rice stubble. Therefore, backscatter coefficient characteristics vary between double-season rice and RR. Additionally, the ratoon crop persistently thrives from rice stubble, with a growth cycle shorter than that of second-season rice in double-cropping systems. In summary, the distinguishability between RR and other paddy rice types contributes to the enhancement of the classification performance of the proposed model.

This study directly compared the performance of proposed model with the performance of other previously reported models for RR mapping (Table 2). It is feasible to directly compare different types of methods, phenological information required, data sources, and spatial resolution. However, comparing the performance of different models under the same conditions of location, time, and sample points is more meaningful. Future research will involve in-depth analysis of existing studies and improvement of comparative experiments.

5.2. Advantages and Limitations of RR Mapping Based on SAR Data

Notably, unlike machine learning methods, which require a large number of training samples [41], our proposed method requires only a small number of field survey samples and a phenological calendar to determine the relevant thresholds, enabling accurate extraction of RR. SAR data are not affected by clouds or fog, making its preprocessing relatively straightforward, typically involving only speckle noise reduction and filtering. When optical data serve as the data source for RR mapping, it is necessary to utilize cloud removal interpolation methods and spatiotemporal fusion models, increasing the intricacy of the processing workflow [15]. In this study, the threshold obtained based on the training samples was essentially consistent with the phenology calendar of RR, indicating the strong sensitivity and good performance of SAR data in identifying critical phenological periods. In contrast, images corresponding to the key phenological stages of RR in optical datasets are mostly obtained through spatiotemporal fusion models; thus, the thresholds determined may differ from the RR phenology calendar [15]. However, compared with some previous studies [15], the RR mapping results based on SAR data exhibit a more pronounced “salt and pepper” effect, and there are substantial differences in some regions when compared with mapping results based on optical data.

Therefore, we recommend that in regions characterized by predominantly clear atmospheric conditions, ORS data should be used [42]; however, in cloudy or rainy weather conditions, if ORS imagery is insufficient to construct a complete time series, the utilization of SAR data is advised. In recent years, the blending of ORS and SAR information for the monitoring of paddy rice fields has emerged as a popular research trend, leveraging the complementary strengths of both modalities. This integration combines spectral characteristics, backscatter coefficient features, and textural features, demonstrating superior performance in paddy rice monitoring compared with single-source data inputs [40]. Future research may explore the incorporation of both optical and SAR data for RR mapping.

5.3. Applicability and Directions for Model Improvement

The proposed threshold model is applicable in locations where there is an available RR phenology calendar to determine the critical characteristics of “double peaks” and “double troughs”. Herein, we identified multiple phenological thresholds based on the RR phenology calendar of Yongchuan District. A phenology calendar may not be available in every region, constraining the widespread application of the threshold model. The long-term stability of the phenological calendar faces challenges from variations in meteorological factors and choices made by farmers [43,44,45]. Therefore, a method that does not rely on phenological information would advance the large-scale and rapid mapping of RR in the future.

A shorter revisit cycle of SAR may help to determine more accurate thresholds. Here, we relied on Sentinel-1A data to map the extent of RR in the study region, because Sentinel-1B data were unavailable, which introduced some uncertainty into our RR mapping. The proposed threshold model utilizes Sentinel-1 data as its basis. Nonetheless, this model could theoretically be adapted to other SAR data, including higher-resolution RADARSAT-2 and Gaofen-3 data. Validating the adaptability of this model to alternative SAR datasets remains a potential avenue for future research.

6. Conclusions

SAR data possess the advantage of penetrating clouds and fog, exhibiting immense potential for application in regions prone to cloudy and rainy conditions. This study utilized SAR data from Sentinel-1A images to acquire the characteristics of VH backscatter coefficients during the growth stages of RR. By integrating the phenological calendar, a threshold model for VH backscatter coefficients was constructed, enabling remote sensing monitoring of RR in Yongchuan District, Chongqing. Sentinel-1A data underwent preprocessing using Lee sigma and S-G filtering. The VH backscatter coefficients of RR, single-season rice, water body, urban land, and forests were analyzed to establish the RR threshold model. The thresholds of critical stages were established by analyzing the VH backscatter time-series curve and referring to the phenological calendar of RR. The robustness of the proposed threshold model was validated in Yongchuan District, Chongqing. Our findings indicated an RR planting area of 199.79 km² in the study region, which shows a high level of consistency with agricultural statistics. Additionally, RR was recognized with high accuracy, in which the OA, F1 score, and Kappa coefficient were 90.24%, 0.92 and 0.80. The results suggest that the proposed threshold model has evident advantages in mapping RR in cloudy and foggy regions. It is more scientific to compare the performance of different models under the same experimental conditions. Future research will carry out in-depth comparative experiments to make up for the limitations of this study. Furthermore, future enhancements to this model should prioritize the integration of optical and SAR data without relying on phenological information.