Skip to content
BY 4.0 license Open Access Published by De Gruyter Open Access May 4, 2022

Mapping paddy rice and rice phenology with Sentinel-1 SAR time series using a unified dynamic programming framework

Mo Wang, Jing Wang, Li Chen and Zhigang Du
From the journal Open Geosciences

Abstract

Monitoring rice planting areas and their phenological phases is crucial for yield estimation and informed decision-making. This study proposed a unified method for mapping rice field and rice phenology with a dynamic time wrapping (DTW) distance-based classifier and its variant sub-DTW algorithm using Sentinel-1’s synthetic aperture radar (SAR) VH band. Field samplings were conducted for broad landcover types in one of the areas of interest (AOIs). We implemented a pixel-wise k-nearest neighbor classification model with DTW distance to identify paddy rice pixels. Standard rice phenological profiles of the SAR VH band were defined by ground monitoring of a sample rice field. Based on rice planting maps and the standard phenological profiles, rice phenological phases were estimated by pattern matching strategy with the sub-DTW algorithm. Experiments on six counties in Northeast China presented promising results. The overall producer and user accuracy reached 92.9 and 91.9% for rice mapping, respectively. The mean root mean square error (RMSE) for phenology estimation was 3.5 days. Rice planting and rice phenology maps were generated for the six AOIs. The phenological variances of the AOIs implied the effects of climate and rice cultivars on phenological development.

1 Introduction

Rice is an important staple food for many countries and contributes 19% of the daily human energy supply, ranked first in the cereal category [1]. Efficient and timely mapping of rice cultivation and rice phenology helps sustainable and precision farming of paddy rice, hence benefiting food security. Besides, paddy rice agriculture plays a considerable role in water resource use, climate change, and even disease transmission [2]. Particularly, information on phenological development is a key to crop monitoring since it depicts the actual states of the crops. Therefore, monitoring paddy rice fields and their growth status has drawn extensive research interest.

Phenological phases of paddy rice are generally divided into three phases: vegetative phase, reproductive phase, and ripening phase, which are further classified into ten growth stages [3]. Paddy rice’s phenological development coincide with changes in underlying water and soil conditions and canopies’ biophysical and biochemical attributes, reflecting changes in time-series remote sensing data.

Remote sensing has provided data sources for rice mapping and phenology detection with frequent acquisition and low cost of data for decades. Various remote-sensing data and algorithms have been applied in related publications, either for rice field mapping or for phenology detection. Widely used remote sensing data sources include optical (e.g., Landsat series [4,5], MODIS [6,7], SPOT [8], Sentinel-2 [9], and HJ-1 [10,11,12]) and microwave (e.g., RADARSAT [13], ALOS/PALSAR [14], TerraSAR-X [15,16], and Sentinel-1 [17,18,19]). A rice crop canopy’s reflectance spectrum results from a complex relationship between its biophysical and biochemical properties [20]. Paddy rice is the only crop that grows in wetland conditions. Some unique time-series features can be explored through its growing cycle. Particularly, flooding signals of the early growing stages are often used to identify paddy rice with optical or microwave remote sensing data, e.g., Yin et al. [21], Xiao et al. [22], and Qiu et al. [23]. While using optical remote sensing data, band combinations or vegetation/water index from multi-spectral bands are key indicators input to classification models. With microwave remote sensing, synthetic aperture radar (SAR) data were also widely used for mapping paddy through detecting key phenological phases, e.g., Torbick et al. [14], Zhang et al. [24], Bazzi et al. [25], and Onojeghuo et al. [26].

Phenology estimation is more challenging than mapping paddy rice with satellite remote sensing data due to gaps in time-series remote sensing data. Analogous to general plant phenology estimation, paddy rice phenology detection via time-series remote sensing data typically consists of three major steps: (1) data cleaning and flagging, (2) data smoothing and time-series data reconstruction, and (3) extraction of phenological metrics based on the reconstructed time-series data [27]. Like paddy rice mapping, band combinations of multi-spectral data are reported more sensitive to indicate crop status [28]. Various vegetation indices were adopted for rice phenology characterization, e.g., NDVI [12,29], and combined with LSWI [30], EVI [31], and EVI2 [11]. Besides, physical-based indexes, e.g., LAI [32,33], were used to characterize plant phenology. More broadly, SAR data with different polarizations or combinations were used for phenology retrieval, owing to microwave remote sensing data integrity (all-weather capability). Different polarization schemes from time-series data of TerraSAR-X [15,34], RADARSAT-2 [35,36], and Sentinel-1 [37,38] were employed for rice phenology monitoring. Multi-polarizations and multi-angular data of SAR signals will provide richer information for rice growth monitoring. However, for the Sentinel-1 C-band data source, several studies have suggested that VH polarization is more sensitive to rice growth than VV since VV polarization is more influenced by standing water in fields and the signal is attenuated by the vertical structure of rice [17,39].

Most of the rice phenology monitoring studies with optical or radar remote sensing data treated phenology detection as a classification or curve-fitting problem, which was widely adopted by many rice mapping studies as well [15,32,36]. Some successful rice mapping studies used a phenology-based method to identify rice fields by identifying some key phenological phases of rice, e.g., Dong et al. [40] and Qiu et al. [23]. This methodology implies that a unified methodology framework could potentially be applied for rice mapping and its phenology monitoring simultaneously. However, to the best of our knowledge, no literature has been found to achieve rice mapping and phenology monitoring under the same methodology framework.

In addition, there are some obstacles while using time-series remote sensing data for regional crop monitoring. (1) One’s study area can be large enough to cover several swaths of multiple orbit tracks, which means sensing dates of images that the study area covers are incoherent. This inconsistency could impede the training of a unified time-series model for the whole study area. (2) For crop mapping or monitoring applications, ground sampling is expensive due to the large extent of the study area. Therefore, it is desirable to develop paddy rice monitoring methods that depend on fewer ground samples. We propose a paddy rice mapping and phenology monitoring framework with a unified dynamic programming strategy to tackle these issues. A k-nearest neighbor (k-NN) classification model based on dynamic time wrapping (DTW) distance was first applied to identify rice field pixels from time-series SAR images. Rice phenology was then estimated with a sub-DTW algorithm from a standard rice phenology profile.

2 Materials and methods

2.1 Study area

The study area consists of six sites in Heilongjiang province, the northmost region of China. The province is the biggest high-quality japonica rice producer in China, with a rough rice-growing area of 400 million hectares. The six sites were selected if there was a National Agricultural Meteorological Station (AMS), which records vital rice phenological phases for surrounding rice fields. The records of AMSs provide field survey data and enable the validation of phenology detection from remote-sensing means.

A humid continental climate predominates in Heilongjiang province. Its winter is chill and dry, with an average of −31 to −15°C in January. Summer is warm and rainy, with an average of 18–23°C in July. The annual average rainfall is 400–700 mm, concentrated heavily in summer. Mudanjiang River basin and Sanjiang Plain of the province have ideal climate and hydrology conditions for rice growing. The annual valid accumulated temperature ≥10°C reaches 2,500 degree days with abundant rainfall and 130–150 frost-free days for these rice-growing regions. Single rice cropping is practiced in the study area. Specifically, transplanted rice predominates in the province, while a small quantity of direct seeding has emerged in recent years. Rice-growing season is from April to October. Rice planting starts in early April, then seedlings are transplanted in early May, and grains are harvested in October.

The selected six counties are Wuchang, Ning’an, Fangzheng, Tangyuan, Baoqing, and Hulin (Figure 1). Each county has an AMS that records the rice phenological phase of the region. Within each county, an area of interest (AOI) was created with a 10 km × 10 km square buffer around the AMS point. As a result, six 100 km2 AOIs were created for this study.

Figure 1 
                  Map of the study area in Heilongjiang province. The six counties are shown with red boundaries. RGB satellite images of the six AOIs are shown in the left and right of the figure.

Figure 1

Map of the study area in Heilongjiang province. The six counties are shown with red boundaries. RGB satellite images of the six AOIs are shown in the left and right of the figure.

2.2 Data

2.2.1 SAR satellite data

In this study, we used time-series SAR satellite images for paddy rice mapping and phenology monitoring. European Space Agency provides a free SAR data source with varied polarization modes by the Sentinel-1 constellation. We acquired the Sentinel-1 Level-1 Ground Range Detected (GRD) data product of the six AOIs for 2019. Specifically, merely the cross-polarization VH band was derived from the GRD product to construct a time-series dataset, since VH polarization is more suitable than VV polarization for rice growth monitoring.

The data preprocessing steps of clipping images to the study area and the time frame were finished on Google Earth Engine (GEE) platform. The SAR GRD data itself were undergone several preprocessing steps before distributing on GEE, including:

  1. Apply orbit file,

  2. GRD border noise removal,

  3. Thermal noise removal,

  4. Radiometric calibration, and

  5. Terrain correction (orthorectification).

Through these steps, original GRD data were processed to backscatter coefficient (σ°) in decibels (dB) suitable for input into models. The Sentinel-1 two-satellite constellation offers a 6-day repeat cycle. However, it should be noted that only Sentinel-1 B data were available online, resulting in a 12-day revisit time if without overlapping of orbits in the AOIs. Acquired time series of Sentinel-1 SAR GRD data for each AOI are listed in Table 1.

Table 1

Sentinel-1 SAR GRD images for each AOI

AOI name SAR time series Interval Image number Bands Resolution
Wuchang 2019-01-08 to 2019-12-22 12 days 28 VH 10 m
Ning’an1 2019-01-11 to 2019-12-25 61
2019-01-03 to 2019-12-29
Fangzheng2 2019-01-08 to 2019-12-22 59
2019-01-03 to 2019-12-29
Tangyuan 2019-01-03 to 2019-12-29 31
Baoqing 2019-01-10 to 2019-12-24 32
Hulin 2019-01-10 to 2019-12-24 30

1,2These two AOIs are located within an overlap of two Sentinel-1 B orbits, therefore having two time series of SAR images.

2.2.2 Landcover sampling data

Field sampling was conducted for five broad landcover types, namely, rice fields, constructions (buildings, roads, bare surfaces, and paved surfaces), water bodies (rivers, lakes, and reservoirs), vegetation (forest, shrubs, orchards, and other planted area), and other croplands. Based on field survey photos and GPS positioning data, sample polygons of the five landcover types were delineated using high-resolution satellite images (CNES/Airbus) on the GEE platform. We collected 147 polygons from Wuchang AOI among which 54 were rice fields, 26 were constructions, 24 were waterbodies, 20 were vegetation, and 23 were other croplands. Another 46 rice field polygons were collected from the other 5 AOIs for test purposes.

The collected sample polygons were stored in a vector format and then rasterized into a GeoTIFF format with 10 m × 10 m pixel spacing to match the resolution of the SAR images. The resulting pixel samples for each landcover type are 5,754 for rice fields, 2,815 for construction, 2,602 for waterbodies, 2,772 for vegetation, and 2,501 for other crops.

Pixel samples were split into the training set and the test set. We deliberately selected 1,000 pixels for each land cover type as the training set from merely one AOI, Wuchang, to test out the robustness of the classification model. Another 1,000 pixel samples of rice field from each AOI were selected to form a test set.

2.2.3 AMS observational data

We collected the AMS observational data that logged the Start-of-Season (SOS) for 11 growth stages, namely (1) seeding, (2) germination, (3) three-leaf, (4) transplanting, (5) reviving, (6) tillering, (7) elongation, (8) booting, (9) heading, (10) milk-ripe, and (11) maturation. An example of the AMS data is shown in Table 2. The AMS observation was conducted on a selected area representing the local rice planting system and dominating rice cultivar near the AMS station.

Table 2

Rice phenology log of the AMS in Wuchang

Phenological phase SOS date (YYYY-MM-DD)
Seeding 2019-04-16
Germination 2019-04-24
Three-Leaf 2019-05-06
Transplanting 2019-05-26
Reviving 2019-06-06
Tillering 2019-06-24
Elongation 2019-07-08
Booting 2019-07-24
Heading 2019-08-01
Milk-ripe 2019-08-16
Maturation 2019-09-16

2.2.4 Phenology observations of sample rice field

We selected a sample rice field to conduct continuous phenology observation. Figure 2 illustrates the location of the sample rice field in Wuchang county. Compared with the 11 stages of AMS observations, we focus on 7 vital phenological phases, i.e., transplanting (AMS stages 1–5), tillering, elongation, booting, heading, milk-ripe, and maturation.

Figure 2 
                     Location of the sample rice field for phenology observation.

Figure 2

Location of the sample rice field for phenology observation.

Following the observation standard, observation was repeated every other day since seeding and until rice’s maturation. SOS of stages of germination and leaf development and flowering maturation was determined if 50% of rice plants reached the phase in the field. For stages of tillering, elongation, booting, and heading, the determining ratio was 10%. The SOS of the six stages was recorded as shown in Table 3.

Table 3

Rice phenology log of the sample rice field

Phenological phase SOS date (YYYY-MM-DD)
Transplanting 2019-05-16
Reviving 2019-06-03
Tillering 2019-06-20
Elongation 2019-07-02
Booting 2019-07-21
Heading 2019-07-28
Milk-ripe 2019-08-11
Maturation 2019-09-08

The recorded SOS of phenological phases and the sample field’s mean value of SAR VH backscatter were retrieved to construct a series of standard phenology profiles (Figure 3), i.e., the curves between each pair of neighboring phenological phases. Standard phenology profiles were extracted for the following pattern matching process.

Figure 3 
                     Sample field’s SAR backscatter profile of pixels’ mean VH band value and the corresponding SOS of phenological phases.

Figure 3

Sample field’s SAR backscatter profile of pixels’ mean VH band value and the corresponding SOS of phenological phases.

2.3 Methods

We propose a method for mapping rice planting and estimating phenological phases under a unified framework using the Sentinel-1 VH polarization band (Figure 4). Pixel-wise rice mapping was first conducted before phenology estimation. We utilized rice’s phenological characteristics of radar signals to design a DTW-based k-NN classifier for rice mapping with a broad landcover classification system. The classification was implemented on the VH band. Following several data processing steps, e.g., interpolation to fill date gaps and data smoothing, phenology estimation of SOS was conducted on rice pixels. We applied a sub-DTW algorithm to the time-series SAR VH band of rice field pixels against standard phenology profiles of each stage to match the SOS dates. The estimated SOS dates were then validated with the phenology logs of AMS.

Figure 4 
                  The workflow of proposed method for mapping rice fields and estimating phenological phases. The left part of the chart shows the rice mapping process with a DTW-based k-NN classifier. The right part shows the phenology estimation process with a pattern-matching strategy using the sub-DTW algorithm. The four data sources are also shown in shaded rounded rectangles.

Figure 4

The workflow of proposed method for mapping rice fields and estimating phenological phases. The left part of the chart shows the rice mapping process with a DTW-based k-NN classifier. The right part shows the phenology estimation process with a pattern-matching strategy using the sub-DTW algorithm. The four data sources are also shown in shaded rounded rectangles.

Figure 5 
                  An illustration of DTW on sequences X and Y with different lengths from Müller [44]. In our case, one sequence is the time series of a sample pixel, and another one is the time series of a pixel of the study area to be classified.

Figure 5

An illustration of DTW on sequences X and Y with different lengths from Müller [44]. In our case, one sequence is the time series of a sample pixel, and another one is the time series of a pixel of the study area to be classified.

2.3.1 Rice mapping with DTW-based k-NN classifier

k-NN algorithm is a widely used non-parametric classification method based on a simple majority vote idea. An object is classified by a majority vote of its neighbors, with the object being assigned to the class having the highest frequency of occurrence among its k-NN. This classification method was also applied in other rice mappings [41] and landcover classification [42] studies. In ref. [41], the authors compared several machine learning or deep learning models, including k-NN, for rice field detection with the Sentinel-1 VH band. However, that study only exploited a single band of SAR data and applied Euclidean distance, neglecting that planting and harvesting schedules may vary for farmers. The potential of k-NN in rice mapping could be promoted with additional polarizations and more elaborate distance measurement involved.

Pixel-wise classification for rice mapping only employed time-series SAR images of the rice-growing season (April to October). We first randomly selected 1,000 pixels from each type of landcover, creating a sample set with 5,000 pixels. A pixel-wise k-NN classifier was then run with DTW distance to each pixel of the sample set. The nearest 200 samples were set (k = 200) to determine the prediction. The mathematical details of DTW are explained.

The DTW algorithm uses dynamic programming techniques to find the optimal alignment and the minimal distance between two temporal sequences. The algorithm was initially developed for speech recognition to cope with different speaking speeds [43] (Figure 5).

Given two sequences X = ( x 1 , x 2 , , x N ) of length N and Y = ( y 1 , y 2 , , y M ) of length M, the objective of DTW is to temporally align these two sequences in some optimal sense under certain constraints.

A warping path W maps the elements of X and Y to minimize the distance between them. W is a sequence of grid points (i,j) in the N-by-M grid. The goal is to find an optimal warping path between X and Y, which is defined to be a warping path P = { ( p 1 , q 1 ) , ( p 2 , q 2 ) , , ( p k , q k ) } that has minimal total cost among all possible warping paths (equation (1)). DTW ( X , Y ) distance between the two sequences, X of length N and Y of length M, is the total cost of an optimal (N, M)-warping path P (equation (2)).

(1) D min ( i k , j k ) = min k 1 , j k 1 D min ( i k 1 , j k 1 ) + d ( i k , j k | i k 1 , j k 1 ) ,

(2) DTW ( X , Y ) =   k d ( i k , j k ) .

Taking x = [3, 1, 2, 2, 1], y = [2, 0, 0, 3, 3, 1, 0] as an example, the accumulated cost matrix and warping path can be illustrated in Figure 6. The warping path follows the path with the lowest cost on the grid.

Figure 6 
                     Accumulated cost matrix and warping path (the yellow dash line) for sequence X and sequence Y.

Figure 6

Accumulated cost matrix and warping path (the yellow dash line) for sequence X and sequence Y.

The warping path is determined using a dynamic programming method to align two sequences. The number of possible paths through the grid will be huge. Therefore, it is critical for efficiency to restrict the number of potential warping paths, and so the following constraints arise.

Boundary condition: this restriction guarantees that the warping path starts at the start points of both sequences and ends at their respective endpoints.

(3) i 1 = 1 , i k = n and j 1 = 1 , j k = m .

Monotonicity condition: this condition restricts the time order of points to ensure that alignment will not go back in time.

(4) i t 1 i t  and   j t 1 j t .

Step-size condition: this condition ensures the path transitions to neighboring points (not jumping in time).

(5) i t i t 1 1   and   j t j t 1 1 .

2.3.2 Time-series reconstruction

Sentinel-1 SAR data source has an interval of 6 or 12 days in this study. To achieve the day-level prediction, we first interpolated the original time series to a per-day sequence and then applied smoothing techniques to achieve time-series reconstruction. The first step was implementing a linear interpolation on the time-series SAR VH band for each pixel to produce a daily time series. We then applied data smoothing on the interpolated time series to acquire continuous and smooth dynamics of the SAR data. Many data smoothing techniques, generally decomposition techniques or curve fitting methods, can be potentially applied for time series of remote-sensing observations [45]. Nonetheless, it is suggested that Savitzky–Golay (SG) filter and wavelet transform (WT) filter had a better curve fitting and data smoothing output. We tested the SG filter and WT filter to reconstruct daily time-series SAR data. The filter with the superior result will be selected.

2.3.2.1 SG filter

The SG filter is a low-pass filtering approach that uses a local polynomial regression model to smooth time-series data. The critical process is a convolution by fitting consecutive subsets of adjacent data points with a low-degree polynomial by the linear least-squares function. As SG is sensitive to smoothing window size, we set the window sizes to 13, considering original data gaps (6 or 12 days).

2.3.2.2 WT filter

WT is a signal processing method that decomposes a function into a set of wavelets. This data filtering method was developed to overcome a significant disadvantage of the Fourier transform in that it captures global frequency information.

Inspired by Sakamoto et al. [46], we designed a synthetic input signal of multi-years with 2,048 daily elements and a cycle of 365 days that repeats the time series of 2019. The year’s VH data were used recursively to fill the input array according to the day of the year (DoY). The missing data of the gaps were linearly interpolated. We used the 365 elements of 1,095th to 1,470th of the output array as the reconstructed SAR VH time profile of 2019. In our case, the growing season of the rice in the study area ranges from 110 to 180 days. Three popular types of mother wavelets were tested: Daubechies (order 2–24), Coiflet (1–5), and Symlet (4–15).

2.3.3 Pattern matching with the sub-DTW algorithm

We focused on estimating phenological phases by ascertaining their SOS with a pattern-matching approach. First, we plotted the standard time profile of each phenological phase’s SAR VH signals with the ground observation on the sample field. The standard time profiles were then aligned with the reconstructed time series of SAR VH data using the sub-DTW algorithm. The SOS of a phenological phase was judged as the date where the alignment starts on the reconstructed time series.

The sub-DTW algorithm is a subsequence variant of the original DTW algorithm. The DTW algorithm fixed the alignment’s beginning match as the first elements of sequences X and Y, respectively, and vice versa for the end match. The sub-DTW algorithm allows for omissions at the beginning and the end Y (longer sequence) in the alignment with X, as illustrated in Figure 7.

Figure 7 
                     An illustration of sub-DTW on sequences X and Y by Müller [44]. In this study, sequence X represents the standard time profile of a phenological phase, and sequence Y represents the time series whose phenological phases are to be estimated.

Figure 7

An illustration of sub-DTW on sequences X and Y by Müller [44]. In this study, sequence X represents the standard time profile of a phenological phase, and sequence Y represents the time series whose phenological phases are to be estimated.

As shown in Figure 8, the SOS of a phenological phase was determined as the beginning date of the matched sequence by the sub-DTW algorithm. For the SOS of the last phenological phase (maturation), it was determined as the ending date of the matched sequenced by the standard profile between milk ripe and maturation in Figure 3.

Figure 8 
                     Phenology SOS detection example with the sub-DTW algorithm.

Figure 8

Phenology SOS detection example with the sub-DTW algorithm.

3 Results

3.1 Rice mapping results

We implemented the DTW-based k-NN classifier on the six AOIs and produced rice planting maps for each AOI (Figure 9). The maps reveal that majority of the agricultural land use was in rice planting for Wuchang county and Fangzheng county. Rice cultivation was less prominent in their agricultural land use for the rest four counties, i.e., Ning’an, Tangyuan, Baoqing, and Hulin.

Figure 9 
                  Rice field maps of the six AOIs by DTW-based k-NN classifier.

Figure 9

Rice field maps of the six AOIs by DTW-based k-NN classifier.

The evaluation of rice mapping was conducted on 1,000 test samples of each AOI with metrics of producer accuracy (PA) and user accuracy (UA). Table 4 shows the evaluation results. Both PA and UA reached a high level for rice mapping in the six AOIs. The best classification performance was observed for Wuchang, Ning’an, and Fangzheng. The reason for Wuchang’s high performance is that the training samples were collected in that AOI. Hence, the classification model was more suited for the region. Rice mappings in Ning’an and Fangzheng had an outstanding performance that can be explained by that they have more time-series images for these two AOIs’ input into the classifier. Therefore, more distinctions between rice fields and other land cover types can be captured by the DTW-based classifier.

Table 4

Evaluation of rice mapping on test samples

AOI Producer accuracy User accuracy Number of images used
Wuchang 0.944 0.938 28
Ning’an 0.933 0.917 61
Fangzheng 0.939 0.925 59
Tangyuan 0.918 0.904 31
Baoqing 0.923 0.907 32
Hulin 0.916 0.920 30
Mean 0.929 0.919

3.2 Timeseries SAR data reconstruction

Following the rice mapping procedure, we implemented the two data smoothing techniques to the time series of identified rice field pixels. The original time series of SAR bands of rice field pixels were reconstructed into daily values throughout 2019.

Figure 10 shows the comparison of original, SG-filtered, and WT-filtered curves. The smoothed curve from the WT filter had a better fitting with the original data points. Some of the inflection points were better restored by the WT filter. Hence, the WT filter was selected for phenology estimation.

Figure 10 
                  Time series of original SAR VH backscatter coefficient in dB and smoothed ones with SG (window size 13) and WT (db6).

Figure 10

Time series of original SAR VH backscatter coefficient in dB and smoothed ones with SG (window size 13) and WT (db6).

3.3 Phenology estimation

Pixel-wise SOS estimation of the eight phenological phases was accomplished with the sub-DTW algorithm. The estimation results were evaluated with AMS loggings. Considering the AMS’s data acquisition standard, the AMS only reflects its nearby rice fields’ phenological phases. Accordingly, validation for our phenology estimation was merely conducted on rice fields within 1 km of the AMS. The mean RMSE was computed for each AOI (Table 5). Phenology estimation achieved relatively high accuracy with RMSE less than 4 days generally. Transplanting, reviving, and heading had relatively higher estimation accuracy due to their prominent inflections in SAR signals. Among the six AOIs, Ning’an and Fangzheng had the least RMSE, which can be explained by more SAR images being involved than other AOIs.

Table 5

Evaluation of phenology estimation within 1 km radius

Phenological phase Mean RMSE of SOS estimation (days)
Wuchang Ning’an Fangzheng Tangyuan Baoqing Hulin
Transplanting 2.4 2.0 2.1 2.6 3.5 3.3
Reviving 2.6 2.3 2.7 3.6 2.7 4.1
Tillering 3.2 2.9 3.5 3.9 2.9 4.2
Elongation 3.5 3.3 3.8 4.7 5.0 5.6
Booting 3.8 3.5 2.6 3.6 4.8 3.9
Heading 2.7 2.6 3.2 3.5 3.9 4.8
Milk-ripe 3.3 3.0 2.9 4.1 4.7 3.8
Maturation 4.1 3.9 3.3 3.8 4.1 4.8
Mean 3.2 2.9 3.0 3.7 3.9 4.3

We chose transplanting, heading, and maturation among the eight phenological phases to examine the spatial variance within and between the AOIs. Estimated transplanting days in DoY for the six AOIs are illustrated as shown in Figure 11. Estimated transplanting days ranged from DoY 122 (May 2, 2019) to DoY 139 (May 19, 2019). Wuchang, Ning’an, and Tangyuan generally had earlier transplanting days than the other three AOIs. This could be explained by Wuchang’s and Ning’an’s lower latitude, which leads to proper air temperature at an earlier date. Tangyuan has a higher latitude but an earlier transplanting day. The possible reason could relate to the particular rice variant. In the meantime, Wuchang and Ning’an showed a more evident variance in transplanting days, indicating a more lenient time window for transplanting. The north parts of these two regions showed earlier transplanting days than the south internally.

Figure 11 
                  Maps of the transplanting day (DoY) for the six AOIs.

Figure 11

Maps of the transplanting day (DoY) for the six AOIs.

Regarding the heading stage, the SOS ranged from DoY 203 (July 22, 2019) to DoY 215 (August 3, 2019), as shown in Figure 12. Wuchang and Ning’an entered the heading stage earlier than the other four AOIs (lighter green than others). A more minor variance could be interpreted between rice fields within each AOI compared to that of transplanting.

Figure 12 
                  SOS maps of the heading stage (DoY) for the six AOIs.

Figure 12

SOS maps of the heading stage (DoY) for the six AOIs.

Figure 13 shows the SOS of the mature stage. The SOS ranged from DoY 250 (September 7, 2019) to DoY 264 (September 21, 2019). The north part of Wuchang and Ning’an showed earlier mature days than the rest regions in general.

Figure 13 
                  SOS maps of the mature stage (DoY) for the six AOIs.

Figure 13

SOS maps of the mature stage (DoY) for the six AOIs.

4 Discussion

The methods for rice mapping or rice phenology estimation with SAR remote sensing data have been widely investigated. It is worth noting that many attempts have achieved satisfying results. In particular, some rice mapping with SAR data utilized phenology traits, for instance, flooding signals during the early growing stage. Such study cases include Clauss et al. [47], Son et al. [48], and Nguyen et al. [17]. Rule-based and machine learning algorithms were prevalently applied methods. We took the same idea of using the unique rice phenology characteristic on radar signals. Instead of utilizing some particular phenological phase, we depicted radar backscattering profiles of the growing season and applied a simple k-NN classification algorithm with DTW distance. The DTW distance resolves the problem of differences in cultivation schedule by time alignment. This metric was also employed by some studies, e.g., Guan et al. [49]. The authors took a threshold on DTW distances of NDVI to discriminate between rice field and non-rice. Our adoption of DTW distance to a more sophisticated classifier would potentially yield a better result.

Rice phenology estimation, on the other hand, was generally treated as the next process on the premise of rice field maps. SAR remote sensing data were also widely used from the earlier TerraSAR-X [34,50,51,52] and RadarSat-2 [36,53,54] to more recent Sentinel-1 [37,55,56]. The authors either relied on certain indexes to detect the maximum/minimal/inflection point on smoothed time series or treated phenology developing progress as a dynamic system and applied Kalman filter or particle filter to estimate the phenology as a system state. These methodologies were proven to be feasible in some of the study cases. However, the first class of methods that detect certain points on curves have its limitation for some mid-term phenology phases, e.g., tillering and booting, in which no evident inflection occurs. The second has its advantage in identifying each phenological phase at any time using prior information from previous observation (remote sensing images). In the meantime, the dynamic system model could bring uncertainties to estimating the SOS of phenological phases. Based on the basic principle that the backscatter signatures of each frequency band have unique patterns that depend on the growth of the rice canopy. With pattern matching of each phenological phase, our approach was straightforward.

Many studies ignored the crop calendar within the study regions. For rice planting, the beginning of a crop cycle is determined by a water distribution scheme. Either rice mapping or phenology estimation shall take this circumstance into account while designing a classification model. The DTW distance could match different crop calendars to avoid this issue, and the pattern-matching strategy was inspired by phenology estimation with the sub-DTW algorithm. Moreover, the pattern matching strategy improves the model versatility for both rice mapping and phenology estimation. Rice cultivation with different crop calendars or different rice varieties would result in uniform or nearly uniform profiles in shape in SAR backscatters. Therefore, in this study, the field samplings for landcover or phenology in one AOI could also apply to the other five AOIs.

The rice mapping and phenology estimation evaluations showed promising performance of the proposed method. The overall PA and UA reached 92.9 and 91.9% for rice mapping, respectively. The mean RMSE for phenology estimation was 3.5 days. Nonetheless, time series that consist of more images had even better performance. This phenomenon was likewise to phenology estimation. Specifically, Ning’an and Fangzheng witnessed the best performance for rice mapping and phenology estimation. Apart from the climate factor, explanations for the phenology differences between the AOIs should furtherly reflect on rice varieties cultivated in different AOIs.

This behavior suggests the importance of filling the data gap in time series. If both Sentinel-1 A and B were available for the AOIs, improvements would take place for rice mapping and phenology estimation. Another noteworthy issue of the proposed method was the neglect of land parcel (rice field) identification. The whole analysis was conducted pixel wisely to test the ability of the method. However, high precise land parcel information would be beneficial to reduce the “salt-and-pepper” effect and furtherly improve phenology estimation in accuracy and speed since land parcels are to be treated as the analytical unit. Super-pixel segmentation with higher-resolution images could be a possible solution to gain land parcel data in the application.

5 Conclusion

This study proposed a unified method for mapping rice field and rice phenology with a DTW distance-based classifier and its variant sub-DTW algorithm. First, field samplings were conducted for broad landcover types in one of the AOIs. We implemented a pixel-wise k-NN classification model with DTW distance to identify paddy rice pixels. Standard rice phenological profiles of the SAR VH band were defined by ground monitoring of a sample rice field. Based on rice planting maps and the standard phenological profiles, rice phenological phases were estimated by pattern matching strategy with the sub-DTW algorithm. Experiments on rice planting regions in Northeast China presented promising results. Rice planting and rice phenology maps were generated for six counties in Heilongjiang province, China. The phenological variances of the AOIs implied the effects of climate and rice cultivars on phenological development.

Acknowledgments

This research was funded by the Central Public-interest Scientific Institution Basal Research Fund of China, grant number Y2021XC17, and Special Project of National Science and Technology Library, grant number 2021XM45.

  1. Funding information: Central Public-interest Scientific Institution Basal Research Fund of China, grant number Y2021XC17, and Special Project of National Science and Technology Library, grant number 2021XM45.

  2. Author contributions: Conceptualization, Mo Wang and Jing Wang; methodology, Mo Wang; software, Mo Wang; validation, Li Chen; resources, Zhigang Du; data curation, Li Chen; writing – original draft preparation, Mo Wang; writing – review and editing, Jing Wang and Zhigang Du; supervision, Jing Wang; project administration, Jing Wang. All authors have read and agreed to the published version of the manuscript.

  3. Conflict of interest: Authors state no conflict of interest.

References

[1] Elert E. Rice by the numbers: a good grain. Nature. 2014;514(7524):S50–1.10.1038/514S50aSearch in Google Scholar PubMed

[2] Dong J, Xiao X. Evolution of regional to global paddy rice mapping methods: a review. ISPRS J Photogramm Remote Sens. 2016;119:214–7.10.1016/j.isprsjprs.2016.05.010Search in Google Scholar

[3] Kuenzer C, Knauer K. Remote sensing of rice crop areas. Int J Remote Sens. 2013;34(6):2101–39.10.1080/01431161.2012.738946Search in Google Scholar

[4] Jin C, Xiao X, Dong J, Qin Y, Wang Z. Mapping paddy rice distribution using multi-temporal Landsat imagery in the Sanjiang Plain, Northeast China. Front Earth Sci. 2016;10(1):49–62.10.1007/s11707-015-0518-3Search in Google Scholar PubMed PubMed Central

[5] Zhang M, Lin H, Wang G, Sun H, Fu J. Mapping paddy rice using a convolutional neural network (CNN) with Landsat 8 datasets in the Dongting Lake Area, China. Remote Sens. 2018;10(11):1840.10.3390/rs10111840Search in Google Scholar

[6] Gumma MK, Nelson A, Thenkabail PS, Singh AN. Mapping rice areas of South Asia using MODIS multitemporal data. J Appl Remote Sens. 2011;5(1):053547.10.1117/1.3619838Search in Google Scholar

[7] Clauss K, Yan H, Kuenzer C. Mapping paddy rice in China in 2002, 2005, 2010 and 2014 with MODIS time series. Remote Sens. 2016;8(5):434.10.3390/rs8050434Search in Google Scholar

[8] Nguyen TTH, De Bie C, Ali A, Smaling E, Chu TH. Mapping the irrigated rice cropping patterns of the Mekong delta, Vietnam, through hyper-temporal SPOT NDVI image analysis. Int J Remote Sens. 2012;33(2):415–34.10.1080/01431161.2010.532826Search in Google Scholar

[9] Cai Y, Lin H, Zhang M. Mapping paddy rice by the object-based random forest method using time series Sentinel-1/Sentinel-2 data. Adv Space Res. 2019;64(11):2233–44.10.1016/j.asr.2019.08.042Search in Google Scholar

[10] Wang J, Huang J, Gao P, Wei C, Mansaray LR. Dynamic mapping of rice growth parameters using HJ-1 CCD time series data. Remote Sens. 2016;8(11):931.10.3390/rs8110931Search in Google Scholar

[11] Wang J, Huang J-F, Wang X-Z, Jin M-T, Zhou Z, Guo Q-Y, et al. Estimation of rice phenology date using integrated HJ-1 CCD and Landsat-8 OLI vegetation indices time-series images. J Zhejiang Univ Sci B. 2015;16(10):832–44.10.1631/jzus.B1500087Search in Google Scholar PubMed PubMed Central

[12] Pan Z, Huang J, Zhou Q, Wang L, Cheng Y, Zhang H, et al. Mapping crop phenology using NDVI time-series derived from HJ-1 A/B data. Int J Appl Earth Observ Geoinf. 2015;34:188–97.10.1016/j.jag.2014.08.011Search in Google Scholar

[13] Park S, Im J, Park S, Yoo C, Han H, Rhee J. Classification and mapping of paddy rice by combining Landsat and SAR time series data. Remote Sens. 2018;10(3):447.10.3390/rs10030447Search in Google Scholar

[14] Torbick N, Salas WA, Hagen S, Xiao X. Monitoring rice agriculture in the Sacramento Valley, USA with multitemporal PALSAR and MODIS imagery. IEEE J Sel Top Appl Earth Observ Remote Sens. 2010;4(2):451–7.10.1109/JSTARS.2010.2091493Search in Google Scholar

[15] Küçük Ç, Taşkın G, Erten E. Paddy-rice phenology classification based on machine-learning methods using multitemporal co-polar X-band SAR images. IEEE J Sel Top Appl Earth Observ Remote Sens. 2016;9(6):2509–19.10.1109/JSTARS.2016.2547843Search in Google Scholar

[16] De Bernardis CG, Vicente-Guijalba F, Martinez-Marin T, Lopez-Sanchez JM. Estimation of key dates and stages in rice crops using dual-polarization SAR time series and a particle filtering approach. IEEE J Sel Top Appl Earth Observ Remote Sens. 2014;8(3):1008–18.10.1109/JSTARS.2014.2372898Search in Google Scholar

[17] Nguyen DB, Gruber A, Wagner W. Mapping rice extent and cropping scheme in the Mekong Delta using Sentinel-1A data. Remote Sens Lett. 2016;7(12):1209–18.10.1080/2150704X.2016.1225172Search in Google Scholar

[18] Wang M, Wang J, Chen L. Mapping paddy rice using weakly supervised long short-term memory network with time series sentinel optical and SAR Images. Agriculture. 2020;10(10):483.10.3390/agriculture10100483Search in Google Scholar

[19] Supriatna R, Shidiq I, Pratama G, Gandharum L, Wibowo A. Spatio-temporal analysis of rice field phenology using Sentinel-1 image in Karawang Regency West Java, Indonesia. Int J. 2019;17(62):101–6.10.21660/2019.62.8782Search in Google Scholar

[20] Yang C-M, Cheng C-H, Chen R-K. Changes in spectral characteristics of rice canopy infested with brown planthopper and leaffolder. Crop Sci. 2007;47(1):329–5.10.2135/cropsci2006.05.0335Search in Google Scholar

[21] Yin Q, Liu M, Cheng J, Ke Y, Chen X. Mapping paddy rice planting area in northeastern china using spatiotemporal data fusion and phenology-based method. Remote Sens. 2019;11(14):1699.10.3390/rs11141699Search in Google Scholar

[22] Xiao X, Boles S, Liu J, Zhuang D, Frolking S, Li C, et al. Mapping paddy rice agriculture in southern China using multi-temporal MODIS images. Remote Sens Environ. 2005;95(4):480–92.10.1016/j.rse.2004.12.009Search in Google Scholar

[23] Qiu B, Li W, Tang Z, Chen C, Qi W. Mapping paddy rice areas based on vegetation phenology and surface moisture conditions. Ecol Indic. 2015;56:79–86.10.1016/j.ecolind.2015.03.039Search in Google Scholar

[24] Zhang Y, Wang C, Wu J, Qi J, Salas WA. Mapping paddy rice with multitemporal ALOS/PALSAR imagery in southeast China. Int J Remote Sens. 2009;30(23):6301–15.10.1080/01431160902842391Search in Google Scholar

[25] Bazzi H, Baghdadi N, El Hajj M, Zribi M, Minh DHT, Ndikumana E, et al. Mapping paddy rice using Sentinel-1 SAR time series in Camargue, France. Remote Sens. 2019;11(7):887.10.3390/rs11070887Search in Google Scholar

[26] Onojeghuo AO, Blackburn GA, Wang Q, Atkinson PM, Kindred D, Miao Y. Mapping paddy rice fields by applying machine learning algorithms to multi-temporal Sentinel-1A and Landsat data. Int J Remote Sens. 2018;39(4):1042–67.10.1080/01431161.2017.1395969Search in Google Scholar

[27] Zeng L, Wardlow BD, Xiang D, Hu S, Li D. A review of vegetation phenological metrics extraction using time-series, multispectral satellite data. Remote Sens Environ. 2020;237:111511.10.1016/j.rse.2019.111511Search in Google Scholar

[28] Galford GL, Mustard JF, Melillo J, Gendrin A, Cerri CC, Cerri CE. Wavelet analysis of MODIS time series to detect expansion and intensification of row-crop agriculture in Brazil. Remote Sens Environ. 2008;112(2):576–87.10.1016/j.rse.2007.05.017Search in Google Scholar

[29] Lin W, Zhang F-C, Jing Y-S, Jiang X-D, Yang S-B, Han X-M. Multi-temporal detection of rice phenological stages using canopy spectrum. Rice Sci. 2014;21(2):108–15.10.1016/S1672-6308(13)60170-5Search in Google Scholar

[30] Sari DK, Ismullah IH, Sulasdi WN, Harto AB. Detecting rice phenology in paddy fields with complex cropping pattern using time series MODIS data. ITB J Sci. 2010;42(2):91–106.10.5614/itbj.sci.2010.42.2.2Search in Google Scholar

[31] Shihua L, Jingtao X, Ping N, Jing Z, Hongshu W, Jingxian W. Monitoring paddy rice phenology using time series MODIS data over Jiangxi Province. China Int J Agric Biol Eng. 2014;7(6):28–36.Search in Google Scholar

[32] Wang C, Zhang Z, Chen Y, Tao F, Zhang J, Zhang W. Comparing different smoothing methods to detect double-cropping rice phenology based on LAI products–a case study in the Hunan province of China. Int J Remote Sens. 2018;39(19):6405–28.10.1080/01431161.2018.1460504Search in Google Scholar

[33] Boschetti M, Busetto L, Ranghetti L, Haro JG, Campos-Taberner M, Confalonieri R, editors. Testing multi-sensors time series of lai estimates to monitor rice phenology: Preliminary results. IGARSS 2018-2018 IEEE International Geoscience and Remote Sensing Symposium. IEEE; 2018.10.1109/IGARSS.2018.8518494Search in Google Scholar

[34] Lopez-Sanchez JM, Cloude SR, Ballester-Berman JD. Rice phenology monitoring by means of SAR polarimetry at X-band. IEEE Trans Geosci Remote Sens. 2011;50(7):2695–709.10.1109/TGRS.2011.2176740Search in Google Scholar

[35] Lopez-Sanchez JM, Vicente-Guijalba F, Ballester-Berman JD, Cloude SR. Polarimetric response of rice fields at C-band: analysis and phenology retrieval. IEEE Trans Geosci Remote Sens. 2013;52(5):2977–93.10.1109/TGRS.2013.2268319Search in Google Scholar

[36] He Z, Li S, Wang Y, Dai L, Lin S. Monitoring rice phenology based on backscattering characteristics of multi-temporal RADARSAT-2 datasets. Remote Sens. 2018;10(2):340.10.3390/rs10020340Search in Google Scholar

[37] Yang H, Pan B, Li N, Wang W, Zhang J, Zhang X. A systematic method for spatio-temporal phenology estimation of paddy rice using time series Sentinel-1 images. Remote Sens Environ. 2021;259:112394.10.1016/j.rse.2021.112394Search in Google Scholar

[38] Wali E, Tasumi M, Moriyama M. Combination of linear regression lines to understand the response of Sentinel-1 dual polarization SAR data with crop phenology – case study in Miyazaki, Japan. Remote Sens. 2020;12(1):189.10.3390/rs12010189Search in Google Scholar

[39] Lasko K, Vadrevu KP, Tran VT, Justice C. Mapping double and single crop paddy rice with Sentinel-1A at varying spatial scales and polarizations in Hanoi, Vietnam. IEEE J Sel Top Appl Earth Observ Remote Sens. 2018;11(2):498–512.10.1109/JSTARS.2017.2784784Search in Google Scholar PubMed PubMed Central

[40] Dong J, Xiao X, Menarguez MA, Zhang G, Qin Y, Thau D, et al. Mapping paddy rice planting area in northeastern Asia with Landsat 8 images, phenology-based algorithm and Google Earth Engine. Remote Sens Environ. 2016;185:142–54.10.1016/j.rse.2016.02.016Search in Google Scholar PubMed PubMed Central

[41] Crisóstomo de Castro Filho H, Abílio de Carvalho Júnior O, Ferreira de Carvalho OL, Pozzobon de Bem P, dos Santos de Moura R, Olino de Albuquerque A, et al. Rice crop detection using LSTM, Bi-LSTM, and machine learning models from sentinel-1 time series. Remote Sens. 2020;12(16):2655.10.3390/rs12162655Search in Google Scholar

[42] Thanh Noi P, Kappas M. Comparison of random forest, k-nearest neighbor, and support vector machine classifiers for land cover classification using Sentinel-2 imagery. Sensors. 2018;18(1):18.10.3390/s18010018Search in Google Scholar PubMed PubMed Central

[43] Sakoe H, Chiba S. Dynamic programming algorithm optimization for spoken word recognition. IEEE Trans Acoust Speech Signal Process. 1978;26(1):43–9.10.1016/B978-0-08-051584-7.50016-4Search in Google Scholar

[44] Müller M. Fundamentals of music processing: audio, analysis, algorithms, applications. Switzerland: Springer; 2015.10.1007/978-3-319-21945-5Search in Google Scholar

[45] Kandasamy S, Baret F, Verger A, Neveux P, Weiss M. A comparison of methods for smoothing and gap filling time series of remote sensing observations–application to MODIS LAI products. Biogeosciences. 2013;10(6):4055–71.10.5194/bg-10-4055-2013Search in Google Scholar

[46] Sakamoto T, Yokozawa M, Toritani H, Shibayama M, Ishitsuka N, Ohno H. A crop phenology detection method using time-series MODIS data. Remote Sens Environ. 2005;96(3–4):366–74.10.1016/j.rse.2005.03.008Search in Google Scholar

[47] Clauss K, Ottinger M, Künzer C. Mapping rice areas with Sentinel-1 time series and superpixel segmentation. Int J Remote Sens. 2018;39(5):1399–420.10.1080/01431161.2017.1404162Search in Google Scholar

[48] Son N-T, Chen C-F, Chen C-R, Toscano P, Cheng Y-S, Guo H-Y, et al. A phenological object-based approach for rice crop classification using time-series Sentinel-1 synthetic aperture radar (SAR) data in Taiwan. Int J Remote Sens. 2021;42(7):2722–39.10.1080/01431161.2020.1862440Search in Google Scholar

[49] Guan X, Huang C, Liu G, Meng X, Liu Q. Mapping rice cropping systems in Vietnam using an NDVI-based time-series similarity measurement based on DTW distance. Remote Sens. 2016;8(1):19.10.3390/rs8010019Search in Google Scholar

[50] Lopez-Sanchez JM, Ballester-Berman JD, Hajnsek I. First results of rice monitoring practices in Spain by means of time series of TerraSAR-X dual-pol images. IEEE J Sel Top Appl Earth Observ Remote Sens. 2010;4(2):412–22.10.1109/JSTARS.2010.2047634Search in Google Scholar

[51] Yuzugullu O, Erten E, Hajnsek I. Rice growth monitoring by means of X-band co-polar SAR: feature clustering and BBCH scale. IEEE Geosci Remote Sens Lett. 2015;12(6):1218–22.10.1109/LGRS.2015.2388953Search in Google Scholar

[52] Yuzugullu O, Marelli S, Erten E, Sudret B, Hajnsek I. Determining rice growth stage with X-band SAR: a metamodel based inversion. Remote Sens. 2017;9(5):460.10.3390/rs9050460Search in Google Scholar

[53] Yang Z, Shao Y, Li K, Liu Q, Liu L, Brisco B. An improved scheme for rice phenology estimation based on time-series multispectral HJ-1A/B and polarimetric RADARSAT-2 data. Remote Sens Environ. 2017;195:184–201.10.1016/j.rse.2017.04.016Search in Google Scholar

[54] He Z, Li S, Lin S, Dai L, editors. Monitoring rice phenology based on Freeman-Durden decomposition of multi-temporal radarsat-2 data. IGARSS 2018-2018 IEEE International Geoscience and Remote Sensing Symposium. IEEE; 2018.10.1109/IGARSS.2018.8517621Search in Google Scholar

[55] Ramadhani F, Pullanagari R, Kereszturi G, Procter J. Automatic mapping of rice growth stages using the integration of SENTINEL-2, MOD13Q1, and SENTINEL-1. Remote Sens. 2020;12(21):3613.10.3390/rs12213613Search in Google Scholar

[56] Supriatna R, Wibowo A, Shidiq IPA, Pratama GP, Gandharum L. Spatio-temporal analysis of rice field phenology using Sentinel-1 image in Karawang Regency West Java, Indonesia. Int J. 2019;17(62):101–6.10.21660/2019.62.8782Search in Google Scholar

Received: 2021-11-25
Revised: 2022-02-13
Accepted: 2022-04-08
Published Online: 2022-05-04

© 2022 Mo Wang et al., published by De Gruyter

This work is licensed under the Creative Commons Attribution 4.0 International License.