Rice Fields Mapping in Fragmented Area Using Multi-Temporal HJ-1A/B CCD Images

Abstract

Rice is one of the most important crops in the world; meanwhile, the rice field is also an important contributor to greenhouse gas methane emission. Therefore, it is important to get an accurate estimation of rice acreage for both food production and climate change related studies. The eastern plain region is one of the major single-cropped rice (SCR) growing areas in China. Subjected to the topography and intensified human activities, the rice fields are generally fragmented and irregular. How remote sensing can meet this challenge to accurately estimate the acreage of the rice in this region using medium-resolution imagery is the topic of this study. In this study, the applicability of the Chinese HJ-1A/B satellites and a two-band enhanced vegetation index (EVI2) was investigated. Field campaigns were carried out during the rice growing season and ground-truth data were collected for classification accuracy assessments in 2012. A stepwise classification strategy utilizing the EVI2 signatures during key phenology stages, i.e., the transplanting and the vegetative to reproductive transition phases, of the SCR was proposed, and the overall classification accuracy was 91.7%. The influence of the mixed pixel and boundary effects to classification accuracy was also investigated. This work demonstrates that the Chinese HJ-1A/B data are suitable data source to estimating SCR cropping area under complex land cover composition.

1. Introduction

Rice is one of the most important crops in the world and provides the main source of energy for more than half of the world population [1]. Additionally, the seasonally flooded rice fields contribute about 5%–19% of total global methane emission, an important greenhouse gas source, to the atmosphere [2,3]. China produced about one third of the world’s rice on about one fifth of the world’s paddy rice land [4]. During the past two decades, the arable land in China declined at a speed of 0.25 million hectares per year [5]. The trend was more obvious in the eastern plain region of China, where intensified human activities have changed the land use and land cover (LULC) patterns dramatically in the last decades. This region has long been one of the major rice growing areas in China, and the cultivar is dominated by the single-cropped rice (SCR). From the food safety, ecological and policy making points of view, a timely and efficient monitoring and mapping of rice cropping area is critical [6,7]. Conventionally, the local government usually estimates the cropping area of rice by field survey; however, it is time-consuming and costly. As a powerful alternative, remote sensing has proved its effectiveness in estimating rice cropping areas from regional to global scales [8,9,10].

In the literature, many different kinds of optical remote sensing data, e.g., the Advanced Very High Resolution Radiometer (AVHRR), Moderate Resolution Imaging Spectroradiometer (MODIS), SPOT VEGETATION and Landsat-MSS, and techniques have been applied in rice cropping area estimating practices [7,8,11,12,13]. The data mentioned above have demonstrated advantages in rice monitoring at regional to global scales due to wide range of coverage and relative long data archiving. However, coarse resolution satellite data is not suitable for precise rice crop mapping in the eastern plain region of China because the rice fields in this region are relatively small, irregular, and fragmented by well-developed roads and dense water networks, and generally mixed with other land cover types. As a consequence, the mixed-pixel problem is prominent and induces temporal uncertainty in discriminating the spectral signatures of rice and the other land cover types [6].

Middle to high spatial resolution satellite data, e.g., Landsat TM/ETM+/OLI, SPOT and China Brazil Earth Resources Satellite (CBERS), are promising in capturing small patches of crop fields [14,15]. However, the cost and relatively long revisit cycles partially offset their advantages in spatial resolution. Specifically, the cloud cover during monsoon season, which is partially overlapped with the major growing season of the SCR, makes it more difficult to obtain qualified remote sensing imageries [16,17]. For applications where the rice phenology information is critically needed, the satellite data with acceptable spatial resolution and more frequent revisit cycle should be more desirable.

The small sun-synchronous satellites for environment and disaster monitoring and forecasting (HJ-1A/B) of China were launched in 2008. HJ-1A/B satellites have a spatial resolution of 30 m and a revisit cycle of four days (the revisit cycle of the constellation is 2 days), with imaging swath of 700 km. The CCD camera onboard HJ-1A/B includes four bands, i.e., blue, green, red and near-infrared, and the spectral range is 0.43–0.90 μm. HJ-1A/B CCD data have been applied in rice area estimation [18,19] and yield prediction [20]. In this study, however, it is our interest to explore the potential of using HJ-1A/B data to extract small, irregular SCR growing area in the eastern plain region of China, where the mixed-pixel problem is serious as mentioned above. Specifically, it is our interest to take advantage of its high revisit feature of the HJ-1A/B data to capture the key phenology spectral signatures of the SCR to facilitate the classification.

The unique physical feature of rice fields and the phenology of the SCR may provide valuable information for remote sensing classification. The rice grows on flooded soils, and the rice fields are a mixture of rice plant and open water during the transplanting and early period of the growing season [21]. As new leaves and tillers emerged, there is an accelerated increase in canopy height and leaf area of the rice. About 50 to 60 days after transplanting, the rice canopy would cover most of the surface area [22], but the leaf area is still increased till the heading stage. After that, the leaf area of rice starts to decease and the leaf color turns to yellow until the ripening and harvest stages. By using time series remote sensing images, the combined field and phenology features of rice, which differentiate the rice field from the other land cover types, may increase the classification accuracy.

To minimize the interference of external environmental factors, various vegetation indices (VIs) are commonly used in practice [23,24,25]. For example, the well-recognized normalized difference vegetation index (NDVI) [26] has been testified to be closely correlated with leaf area, biomass, percent ground cover and crop productivity [27,28,29,30]. Due to the saturation effect, however, NDVI may fail to capture the difference in well-vegetated areas, compared with the enhanced vegetation index (EVI) [31]. In practice, the time series signatures of NDVI and EVI derived from the MODIS and SPOT data had been used to map the area, species (single, early, and late), and key phenologies of rice [13,18,32]. Recently, a novel VI, i.e., the 2-band EVI (EVI2), has been proposed and testified to be comparable with the traditional EVI, and more importantly, it may achieve greater consistencies across sensors because only 2 bands are involved, as compared with 3 bands in EVI [33,34].

In addition to the data used, it is of critical importance to select appropriate classification method to properly mapping the rice fields. It is our interest to compare the classification efficiencies of the commonly used parametric and nonparametric classification algorithms, i.e., the maximum likelihood classifier (MLC) and support vector machines (SVM), with a two-step classification method proposed in this study and specifically designed to classify the rice fields from the other land cover types. The MLC is one of the most commonly used classification techniques [35,36,37]. It is a parametric classification algorithm with the assumption that the class signatures are normally distributed. The SVM is a nonparametric classifier, which projects the training data in the input space into a high dimensional space using a kernel function where the classes are linearly separable [38]. The SVMs have no limitation about the probability distribution forms of the class signature, but its performance largely depends on the kernel used, the parameter choice for the specific kernel, and the method used to generate the SVM [39,40,41].

Our study aimed to investigate the capability of EVI2 in SCR growth monitoring, and to test the feasibility of using HJ-1A/B CCD data to estimate the SCR growing area in the eastern plain region of China. For this purpose, we proposed a simple but effective classification method, which makes use of the time series HJ-1A/B imageries and the specific signatures of EVI2 (including its 1st derivative) at key phenology stages of the SCR. An extensive field campaign was carried out for verification simultaneously. We compared the effectiveness of this method with the parametric and nonparametric classification algorithms, namely MLC and SVM. We also discussed the influence of the mixed-pixel which was typical in the study area and may affect the classification accuracy.

2. Data and Methods

2.1. Study Area

Deqing County lies in the west of Hangjiahu Plain, with mean annual temperature ranging between 13 °C and 16 °C and annual precipitation of 1379 mm (Figure 1). The plain areas mainly distribute at the eastern Deqing, with the altitudes ranging from 4 m along the Beijing-Hangzhou Grand Canal to 721 m on the Tianmu Mountains. Deqing County is part of the SCR growing region in the water network area of north Zhejiang [42], where countless lakes, ponds and winding rivers scattered throughout this region, with the addition of well-developed road networks, leading to fragmented patches of irregular crop land plots. Deqing has a total area of 936 km², and the SCR area in Deqing accounts for more than 91% of the major crop areas according to the statistical data of local agriculture department. The SCR fields mainly concentrate in the eastern regions of Deqing with average elevation less than 20 m.

Figure 1. HJ-1A CCD false-color composite image of the study area on 19 November 2012.

Figure 1. HJ-1A CCD false-color composite image of the study area on 19 November 2012.

2.2. Field Campaigns

To facilitate the remote sensing classification and verification, a continuous field campaign was carried out to record the phenologies of the SCR. Additionally, five field sites named A to E were also selected at the east of Deqing County. All the sites were larger than 1 km², and were surveyed using a handheld GPS receiver (Trimble Juno-SB). For each land cover patch, the boundary and the corresponding land cover type were recorded. The land cover types were classified as rice, trees, water bodies, economic crops and other nonvegetated areas. The vector format maps of the five field sites in 2012 were shown in Figure 2. These maps were then reclassified into SCR and non-rice area and converted into raster format at 30 m resolution as ground-truth data for accuracy assessment.

Figure 2. Vector maps and geo-locations of the five field sites.

Figure 2. Vector maps and geo-locations of the five field sites.

2.3. Remote Sensing Data

HJ-1A/B data from 17 May to 5 December 2012 over the study area were collected for time-series VI analysis and downloaded from the China Center for Resources Satellite Data and Application. The sensor characteristics are presented in Table 1. Total 14 HJ-1A/B images with cloud cover less than 10% during the key phenology periods of SCR were selected for the following classification procedures (Table 2).

Table 1. Technical specification of HJ-1-A/B CCD and ZY1-02C P/MS sensors.

Satellite	Payload	Band No.	Spectral Range (µm)	Nadir Spatial Resolution (m)	Swath Width (km)	Repetition Cycle (day)
HJ-1A/B	Multispectral CCD camera	1	0.43–0.52	30	360 (700 for two)	4
		2	0.52–0.60	30
		3	0.63–0.69	30
		4	0.76–0.90	30
ZY1-02C	P/MS camera	1	0.51–0.85	5	60	3–5
		2	0.52–0.59	10
		3	0.63–0.69	10
		4	0.77–0.89	10

Table 1. Technical specification of HJ-1-A/B CCD and ZY1-02C P/MS sensors.

Satellite	Payload	Band No.	Spectral Range (µm)	Nadir Spatial Resolution (m)	Swath Width (km)	Repetition Cycle (day)
HJ-1A/B	Multispectral CCD camera	1	0.43–0.52	30	360 (700 for two)	4
		2	0.52–0.60	30
		3	0.63–0.69	30
		4	0.76–0.90	30
ZY1-02C	P/MS camera	1	0.51–0.85	5	60	3–5
		2	0.52–0.59	10
		3	0.63–0.69	10
		4	0.77–0.89	10

Table 2. Dates of the selected HJ-1A/B CCD images, field campaigns, and the corresponding SCR phenology stages.

NO.	Satellite	Date	Field Campaign Date	Phenology Stage
1	HJ-1B	2012/05/17	/	Fallow
2	HJ-1B	2012/05/28	/	Site preparation
3	HJ-1A	2012/06/29	2012/06/29	Sowing-transplanting
4	HJ-1B	2012/07/05	/	Vegetative stage
5	HJ-1B	2012/07/19	2012/07/20	Vegetative stage (tillering)
6	HJ-1A	2012/07/29	2012/07/30	Vegetative stage (maximum tiller number)
7	HJ-1A	2012/08/17	2012/08/15	Reproductive stage (ear differentiation)
8	HJ-1B	2012/09/02	2012/08/31	Reproductive stage (heading)
9	HJ-1B	2012/09/18	2012/09/16	Reproductive stage (panicle initiation and flowering)
10	HJ-1B	2012/09/29	2012/09/25	Reproductive stage
11	HJ-1B	2012/10/10	2012/10/13	Ripening stage (grain filling)
12	HJ-1A	2012/10/23	2012/10/27	Ripening stage (milk)
13	HJ-1B	2012/11/06	/	Ripening stage
14	HJ-1A	2012/11/19	2012/11/18	Harvest/fallow

Table 2. Dates of the selected HJ-1A/B CCD images, field campaigns, and the corresponding SCR phenology stages.

NO.	Satellite	Date	Field Campaign Date	Phenology Stage
1	HJ-1B	2012/05/17	/	Fallow
2	HJ-1B	2012/05/28	/	Site preparation
3	HJ-1A	2012/06/29	2012/06/29	Sowing-transplanting
4	HJ-1B	2012/07/05	/	Vegetative stage
5	HJ-1B	2012/07/19	2012/07/20	Vegetative stage (tillering)
6	HJ-1A	2012/07/29	2012/07/30	Vegetative stage (maximum tiller number)
7	HJ-1A	2012/08/17	2012/08/15	Reproductive stage (ear differentiation)
8	HJ-1B	2012/09/02	2012/08/31	Reproductive stage (heading)
9	HJ-1B	2012/09/18	2012/09/16	Reproductive stage (panicle initiation and flowering)
10	HJ-1B	2012/09/29	2012/09/25	Reproductive stage
11	HJ-1B	2012/10/10	2012/10/13	Ripening stage (grain filling)
12	HJ-1A	2012/10/23	2012/10/27	Ripening stage (milk)
13	HJ-1B	2012/11/06	/	Ripening stage
14	HJ-1A	2012/11/19	2012/11/18	Harvest/fallow

To assist the selection of training samples for classification, the Chinese Resource-1 02C satellite (ZY1-02C), which provides multispectral and panchromatic images at 10 m and 5 m spatial resolutions, respectively, was used as an auxiliary data source (Table 1). The multispectral and panchromatic images of ZY1-02C were fused to facilitate location identity in field campaigns and visual interpretation.

All the HJ 1-A/B and ZY1-02C images were geometrically corrected using the Second National Soil Survey Vector Map (scale 1:10,000), and the Root Mean Square Error (RMS error) was less than one pixel (30 m). Additionally, the radiometric calibration and atmospheric correction of the HJ 1-A/B CCD data were performed, respectively. Figure 3 showed the images of HJ-1A/B CCD and ZY1-02C of field site B at different phenology stages of the SCR.

The remote sensing classification system (five land cover types) was same as the one used in the field campaign. The training samples were randomly located and visually interpreted from the ZY1-02C fused image (5 m in spatial resolution). The class separability of the training data set was analyzed using the Jeffries-Matusita (J-M) distance metric between classes [43,44]. A larger J-M distance indicates more distinct distributions between two classes. The training data were modified if the J-M distance was close to 2 between rice and the other land cover types [45,46]. The final set of training samples were 800 pixels in total. There were 354 training pixels for rice, 92 training pixels for trees, 195 training pixels for water bodies, 76 training pixels for economic crops and 83 training pixels for other nonvegetated areas. The vector data of the 5 field sites were rasterized into 30 m resolution as ground-truth data for accuracy assessment.

Figure 3. HJ-1A/B CCD and ZY1-02C false color images of the field site B at different phenology stages of the SCR: (a) to (c) were acquired from HJ-1A/B on 29 June (sowing-transplanting stage), 29 July (vegetative stage), and 2 September (reproductive stage) 2012, respectively; and (d) ZY1-02C on 19 February 2012.

Figure 3. HJ-1A/B CCD and ZY1-02C false color images of the field site B at different phenology stages of the SCR: (a) to (c) were acquired from HJ-1A/B on 29 June (sowing-transplanting stage), 29 July (vegetative stage), and 2 September (reproductive stage) 2012, respectively; and (d) ZY1-02C on 19 February 2012.

2.4. Classification Methods

2.4.1. Characteristics of EVI2 Time-Series Data during SCR Growing Periods

EVI2 may achieve greater consistencies across sensors because only 2 bands are involved. EVI2 is defined as follows [33]:

$$EVI2=2.5\times \frac{{\rho }_{nir}-{\rho }_{red}}{{\rho }_{nir}+2.4\times {\rho }_{red}+1}$$

(1)

where ρ_nir and ρ_red are estimated surface reflectance values for near-infrared and visible red bands (HJ-1A/B CCD bands 4 and 3, respectively).

Since there are nearly always disturbances in optical remote sensing applications caused by unfavorable atmospheric conditions and sun zenith angle changes in year around and show up as undesirable noise [47,48], noise reduction is necessary before further analysis. In this study, we used the Savitzky-Golay (S-G) filters to smooth the EVI2 time-series data. The S-G filters are suitable to smooth the irregular spacing data points, e.g., the time-series HJ-1A/B CCD data used in this study [49]. The S-G filters apply an iterative weighted moving average filter to time series data, with weighting given as a polynomial of a particular degree [50,51], which can be summarized as:

$$g_i=\frac{{\displaystyle {\sum }_{n=-nL}^{nR}{c}_n{f}_{i+n}}}{n}$$

(2)

where f_i represents original value at data point i; g_i is the smoothed value; n is the width of the moving window to perform filtering; nL and nR corresponding to the left and right edges of the signal component. For a specific uneven time-series data in a moving window, c_n is not a constant but a polynomial fitting function, depending on the user’s preference. The fitting function can be defined as quadratic polynomial for a specific f_i:

$$c_n(t)=c_1+c_{2}t+c_3{t}^2$$

(3)

where t corresponds to the day of year in EVI2 time-series.

The S-G filter was implemented using IDL 8.0 programming language to perform an image-based EVI2 time-series filtering for the HJ-1A/B CCD data from 17 May to 5 December 2012 in the study area. In this way, the time-series EVI2 curves of the five land cover types could be used to identify the most critical stages to distinguish between different land cover types.

2.4.2. Single-Cropped Rice Classification Method

In this study, we proposed a classification method which is based on the assumption that the probability distribution functions (PDFs) of the land cover types follow normal distributions [52,53]. For this purpose, we tested the samples’ probability distributions of the five land cover types using the Quantile-Quantile Plot (Q-Q Plot) [54], and it showed that normal distribution assumption was acceptable.

Using the training data set, the mean (μ) and standard deviation (σ) of each land cover type can be obtained, and then we can define the normal distribution function for each land cover type using these two parameters. To properly differentiate one specific land cover type from the others, it is crucial to minimize the overlaps between the target and the neighboring normal PDFs. For two land cover types L₁ and L₂, assuming L₁~N(μ₁, σ₁²) and L₂~N(μ₂, σ₂²), then the intersection between L₁ and L₂ should be as follows [55]:

$$x=\frac{{\sigma }_1{\mu }_2+{\sigma }_2{\mu }_1}{{\sigma }_1+{\sigma }_2}$$

(4)

where if x is out of [μ − 2σ, μ + 2σ] (hereafter μ(σ) can be μ₁(σ₁) or μ₂(σ₂) either), the two classes can be assumed distinguishable; if x is out of [μ − σ, μ + σ] but within [μ − 2σ, μ + 2σ], the two classes are mildly overlapped; if x is within [μ − σ, μ + σ], then the two classes are seriously overlapped. Generally, the two classes can be thought separable if x is out of [μ − σ, μ + σ].

Instead of using the whole growing period dataset, only key phenology stages images (during which the SCR are most differentiable from the other land cover types as explained later) were investigated in SCR field extraction. We used both of EVI2 and its 1st derivative, calculated by three consecutive images, to minimize the probability of mis-classification. For EVI2, we selected the image on 29 June 2012 (transplanting stage), whilst the spectral characteristic of the SCR was similar to water but not to the other land cover types, especially the trees. In addition, we made use of the quick change rate of EVI2, i.e., the 1st derivative of SCR during the vegetative stages (here we used the image on 29 July 2012) to gather further information to refine the classification results [56,57].

2.4.3. Parametric and Nonparametric Classification Algorithms

The MLC assumes that the class signatures are normally distributed and calculates the probabilities of a given pixel belonging to each class. The pixel is assigned to the class with the highest probability [58]. The SVM classifier is a kernel-based machine learning technique; it separates the classes with a decision surface which maximizes the margin between the classes. The success of the SVM depends on how well the process is trained. In this study, a well-known radial basis function (RBF) kernel was used in the SVM [19,38,41].

We applied the MLC and SVM using the same training samples and parameters for each classifier. The multi-temporal HJ-1 CCD data, i.e., six scenes from 2012/06/29 to 2012/09/02 during SCR transplanting to early reproductive stages, were used. The six reflectance/EVI2 imageries were composited and classified using the MLC and SVM separately, and the results were compared with the method proposed in Section 2.4.2.

2.4.4. Classification Accuracy Assessment

To assess classification accuracy, the ground-truth data (30 m resolution) were used. The ground-truth pixel numbers for field site A to E were 2842, 3078, 3186, 3248, and 2350, respectively. The proposed classification results were compared with the local agricultural statistic data in 2012. The user’s and producer’s accuracies, overall accuracy and Kappa statistic were also used to evaluate the SCR classification accuracy among the proposed method and the traditional methods, i.e., MLC and SVM.

2.5. Influence of the Mixed-Pixel

To analyze the relationship between the land cover structure (or fragmentation) and classification accuracy, we calculated the landscape metrics, i.e., class area (CA), percent of landscape (%LAND), patch density (PD), mean patch size (MPS), area-weighted mean shape index (AWMSI), and mean nearest-neighbor distance (MNN), for each ground-truth site at class level using FRAGSTATS to quantify its structure property [59]. Among the landscape metrics used here, CA is a measure of how much of the landscape is composed of a particular land cover type; %LAND is the percent of each land cover type; PD is the number of patches on a per unit area; MPS is the average area of patches for a certain class; AWMSI measures the area-weighted average patch shape; and MNN measures the mean average nearest distance among patches in a class. The landscape indices, e.g., PD, MPS, AWMSI and MNN, can be used to represent the fragmentation of land cover for a specific field site. The larger the values of PD, AWMSI and MNN, the more fragmented the site was; and vice versa for MPS.

To evaluate the influence of the land cover composition of a specific pixel on the classification accuracy, the vector maps of the five ground-truth sites were further divided into cells of size 30 m × 30 m using the gridlines derived from the HJ CCD images (only the SCR fields were kept and the other land cover types were taken as background). We calculated the area proportion of SCR in each cell, and divided the cells, in which the SCR area were greater than zero, into three grades, i.e., 75%–100%, 50%–75%, and <50%, according to the SCR area proportion. For each grade in per site, the proportion of cells, which were classified as rice field in the HJ CCD images, to the total cell number in that specific grade was calculated. We further calculated the number of misclassified pixels, i.e., the commission and omission errors, and analyzed the corresponding spatial distribution of the misclassified pixels. For a specific pixel, the commission error means that the pixel’s SCR area proportion is less than 50% but is classified as rice field; while the omission error means that the pixel contains more than 50% SCR area but is classified as the other land cover types.

3. Results

3.1. Time-Series EVI2 Characteristics

The temporal dynamics of time-series EVI2 of the SCR processed by S-G filters and the other land cover types calculated from HJ-1A/B images were shown in Figure 4 and Figure 5. The EVI2 of water bodies varied slightly over the growing season of SCR in the range of 0.07–0.14. During the transplanting and early part of the SCR growing period, rice fields were flooded and its spectral signature was similar to that of the water bodies. Not surprisingly, the EVI2 value of rice fields was very close to water bodies but obviously lower than that of trees and economic crops on June 29 (DOY = 181, about 10 days after transplanting in 2012). After transplanting, the EVI2 of SCR increased rapidly and maximized at about 0.6 between the ear differentiation and early heading stages, about 75 days after transplanting. Caused by the etiolation and senescence of the SCR leaves, the EVI2 started to decrease after the heading period till harvest.

The EVI2 values of the other nonvegetated areas, including residential areas, roads and bare land, were similar with less fluctuation but relatively higher compared with water bodies. The trees class had relatively high EVI2 values around 0.30 to 0.47. The economic crops were generally planted during a similar period as the SCR were transplanted, but usually have a longer life cycle and relativly small changing rate of EVI2 compared with SCR, especially during the vegetative stages of SCR. During the transplanting period, the water like spectral characteristic of the SCR made its EVI2 signature a little lower than the economic crops.

Figure 4. Time-series EVI2 data fitted by the S-G filters.

Figure 4. Time-series EVI2 data fitted by the S-G filters.

Figure 5. Time-series EVI2 of the five land use types processed by the S-G filters.

Figure 5. Time-series EVI2 of the five land use types processed by the S-G filters.

3.2. Classification Thresholds

The normal distributions of EVI2 on 29 June 2012 and its first derivative on 29 July 2012 of the five land cover types were shown in Figure 6. During the transplanting stage, the rice fields were flooded and the PDF of EVI2 of the SCR was close to that of the water bodies, and it also mixed with the nonvegetated areas and economic crops classes. Therefore, except the trees, which were mildly separable with respect to the economic crops but highly distinguishable from the other classes, the EVI2 signature of SCR was seriously overlapped with the other three land cover types (Figure 6a). Obviously, it is insufficient to identify the SCR just using images in the transplanting period.

Figure 6. Normal distributions of the five land cover types: (a) EVI2 on 29 June 2012; and (b) the 1st derivative of EVI2 on 29 July 2012.

Figure 6. Normal distributions of the five land cover types: (a) EVI2 on 29 June 2012; and (b) the 1st derivative of EVI2 on 29 July 2012.

During the vegetative stages of the SCR, there was a quick increase of EVI2 due to the formation of additional tillers (Figure 4), while the increase rates of the economic crops and trees were not as steep as SCR. The other two land cover types didn’t show obvious changes during this period. In this situation, the 1st derivative of EVI2 based on the image on 29 July 2012, when the maximum tiller number had arrived, demonstrated that the rice fields could be confidently distinguished from the water bodies and the other nonvegetated areas (Figure 6b). As shown in Figure 6a, the rice fields class was distinguishable from the trees using EVI2 signatures during the transplanting stage, and it is mildly separable from the economic crops class. By using Equation (4), one pixel could be classified as rice field if its EVI2 value on 29 June 2012 equal or less than 0.24, whilst its 1st derivative of EVI2 on 29 July 2012 was equal or greater than 0.010.

3.3. Classification Accuracy Assessment

Based on the coupling thresholds of EVI2 and its 1st derivative of the SCR, during the key phenology stages, i.e., the transplanting and the vegetative to reproductive transition phases, the rice fields in Deqing County was classified using HJ-1A/B data in 2012 (Figure 7). The classified rice fields, with an area of about 94.0 km², mainly concentrated in the eastern plain region of the study area with altitudes around 4 m. According to the statistical data of the local agriculture department in 2012, the total acreage of the SCR was 86.4 km², so the relative classification accuracy was about 91.2%.

Figure 7. Classification result of the SCR fields using HJ-1A/B data in 2012.

Figure 7. Classification result of the SCR fields using HJ-1A/B data in 2012.

Table 3. Classification accuracies and Kappa coefficients of the five ground-truth sites for SCR. The last column corresponded to that if all sites were treated as a whole.

Classification Accuracy	A	B	C	D	E	All
Producer’s accuracy (%)	95.43	94.12	90.72	94.98	95.67	94.35
User’s accuracy (%)	75.1	77.29	69.17	86.31	83.42	78.70
Overall accuracy (%)	92.43	90.68	91.93	94.21	93.40	91.68
Kappa coefficient	0.76	0.75	0.71	0.85	0.83	0.79

Table 3. Classification accuracies and Kappa coefficients of the five ground-truth sites for SCR. The last column corresponded to that if all sites were treated as a whole.

Classification Accuracy	A	B	C	D	E	All
Producer’s accuracy (%)	95.43	94.12	90.72	94.98	95.67	94.35
User’s accuracy (%)	75.1	77.29	69.17	86.31	83.42	78.70
Overall accuracy (%)	92.43	90.68	91.93	94.21	93.40	91.68
Kappa coefficient	0.76	0.75	0.71	0.85	0.83	0.79

We compared the classification accuracies of the 5 ground-truth sites (Table 3). The overall classification accuracy and Kappa coefficient for all the sites were 91.68% and 0.79, respectively. For each site, the producer’s and user’s accuracies of rice, overall classification accuracy and Kappa coefficient are listed in Table 3, in which the site D had the highest user’s accuracy and overall classification accuracy (86.31% and 94.21%, respectively), followed by site E, with user’s accuracy of 83.42% and overall accuracy of 93.40%. All the ground-truth sites had producer’s accuracies higher than 90.72% (site C). The accuracy assessment demonstrated a satisfactory result of the proposed classification method.

3.4. Comparison of Classification Methods

The classification accuracy of the proposed method used in this study outperformed the MLCs and SVMs (Table 4 and Table 5). By using EVI2 instead of the reflectance data, the MLC-EVI2 and SVM-EVI2 improved the classification accuracies to certain extent, but not significant, compared with their counterparts, respectively. The MLCs and SVMs also showed better classification accuracies in site D and site E. While site C had a lower classification result compared with the other four field sites.

Table 4. Classification accuracies and Kappa coefficients of the five ground-truth sites for SCR using MLC and SVM methods. The last column is the corresponding results if all sites were treated as a whole.

Classification Methods	Classification Accuracy	A	B	C	D	E	All
MLC-EVI2	Producer’s accuracy (%)	67.77	70.77	60.68	77.01	70.69	62.41
	User’s accuracy (%)	67.08	63.83	61.43	73.25	73.68	60.58
	Overall accuracy (%)	88.85	85.10	84.03	89.17	87.81	84.88
	Kappa coefficient	0.71	0.68	0.68	0.60	0.74	0.65
MLC-Reflectance	Producer’s accuracy (%)	68.60	66.48	52.97	79.93	70.27	62.55
	User’s accuracy (%)	65.23	61.21	56.15	69.49	70.27	57.52
	Overall accuracy (%)	88.42	83.74	81.54	88.58	86.72	83.83
	Kappa coefficient	0.70	0.63	0.64	0.70	0.72	0.65
SVM-EVI2	Producer’s accuracy (%)	82.23	83.38	60.05	83.76	81.07	73.97
	User’s accuracy (%)	50.90	52.43	50.90	59.14	61.16	54.33
	Overall accuracy (%)	83.46	80.17	76.96	85.37	85.34	81.09
	Kappa coefficient	0.63	0.62	0.57	0.55	0.69	0.61
SVM-Reflectance	Producer’s accuracy (%)	84.30	85.82	60.96	85.40	79.10	74.80
	User’s accuracy (%)	50.00	51.63	51.51	58.40	63.87	52.87
	Overall accuracy (%)	82.97	78.23	76.35	84.78	84.27	80.00
	Kappa coefficient	0.63	0.59	0.56	0.55	0.61	0.57

Table 4. Classification accuracies and Kappa coefficients of the five ground-truth sites for SCR using MLC and SVM methods. The last column is the corresponding results if all sites were treated as a whole.

Classification Methods	Classification Accuracy	A	B	C	D	E	All
MLC-EVI2	Producer’s accuracy (%)	67.77	70.77	60.68	77.01	70.69	62.41
	User’s accuracy (%)	67.08	63.83	61.43	73.25	73.68	60.58
	Overall accuracy (%)	88.85	85.10	84.03	89.17	87.81	84.88
	Kappa coefficient	0.71	0.68	0.68	0.60	0.74	0.65
MLC-Reflectance	Producer’s accuracy (%)	68.60	66.48	52.97	79.93	70.27	62.55
	User’s accuracy (%)	65.23	61.21	56.15	69.49	70.27	57.52
	Overall accuracy (%)	88.42	83.74	81.54	88.58	86.72	83.83
	Kappa coefficient	0.70	0.63	0.64	0.70	0.72	0.65
SVM-EVI2	Producer’s accuracy (%)	82.23	83.38	60.05	83.76	81.07	73.97
	User’s accuracy (%)	50.90	52.43	50.90	59.14	61.16	54.33
	Overall accuracy (%)	83.46	80.17	76.96	85.37	85.34	81.09
	Kappa coefficient	0.63	0.62	0.57	0.55	0.69	0.61
SVM-Reflectance	Producer’s accuracy (%)	84.30	85.82	60.96	85.40	79.10	74.80
	User’s accuracy (%)	50.00	51.63	51.51	58.40	63.87	52.87
	Overall accuracy (%)	82.97	78.23	76.35	84.78	84.27	80.00
	Kappa coefficient	0.63	0.59	0.56	0.55	0.61	0.57

Table 5. Landscape indices of the five land cover types in the ground-truth sites.

Site	Class	CA (ha)	%LAND	PD	MPS (ha)	AWMSI	MNN (m)
A	Rice	50.08	25.78	1.56	0.64	1.73	27.41
	Trees	44.99	23.16	2.20	0.45	2.83	19.45
	Water bodies	69.12	35.58	1.22	0.82	2.54	18.03
	Economic crops	2.32	1.19	6.89	0.15	1.40	126.17
	Other nonvegetated areas	27.77	14.29	4.83	0.21	7.76	30.91
B	Rice	60.85	30.25	1.25	0.80	2.11	21.90
	Trees	56.73	28.20	1.73	0.58	3.23	21.38
	Water bodies	48.12	23.92	1.41	0.71	3.86	20.32
	Economic crops	3.75	1.86	5.33	0.19	5.67	80.80
	Other nonvegetated areas	31.72	15.77	4.48	0.22	4.55	17.64
C	Rice	40.50	19.66	1.31	0.76	1.57	22.14
	Trees	45.47	22.07	1.83	0.55	2.56	15.98
	Water bodies	74.08	35.95	0.90	1.11	2.10	20.15
	Economic crops	7.55	3.67	2.25	0.44	1.34	92.04
	Other nonvegetated areas	38.44	18.66	4.63	0.22	7.34	21.30
D	Rice	70.31	36.03	0.95	1.05	1.25	13.85
	Trees	32.12	16.46	1.96	0.51	2.20	13.24
	Water bodies	62.61	32.08	1.15	0.87	1.94	7.00
	Economic crops	5.95	3.05	15.12	0.07	4.38	21.55
	Other nonvegetated areas	24.15	12.38	8.49	0.12	18.53	13.50
E	Rice	55.62	31.15	1.11	0.90	1.40	13.64
	Trees	41.63	23.32	2.45	0.41	1.97	14.60
	Water bodies	43.74	24.50	0.78	1.29	2.60	22.36
	Economic crops	5.88	3.29	5.44	0.18	1.58	65.06
	Other nonvegetated areas	31.67	17.74	5.27	0.19	8.29	19.73

Table 5. Landscape indices of the five land cover types in the ground-truth sites.

Site	Class	CA (ha)	%LAND	PD	MPS (ha)	AWMSI	MNN (m)
A	Rice	50.08	25.78	1.56	0.64	1.73	27.41
	Trees	44.99	23.16	2.20	0.45	2.83	19.45
	Water bodies	69.12	35.58	1.22	0.82	2.54	18.03
	Economic crops	2.32	1.19	6.89	0.15	1.40	126.17
	Other nonvegetated areas	27.77	14.29	4.83	0.21	7.76	30.91
B	Rice	60.85	30.25	1.25	0.80	2.11	21.90
	Trees	56.73	28.20	1.73	0.58	3.23	21.38
	Water bodies	48.12	23.92	1.41	0.71	3.86	20.32
	Economic crops	3.75	1.86	5.33	0.19	5.67	80.80
	Other nonvegetated areas	31.72	15.77	4.48	0.22	4.55	17.64
C	Rice	40.50	19.66	1.31	0.76	1.57	22.14
	Trees	45.47	22.07	1.83	0.55	2.56	15.98
	Water bodies	74.08	35.95	0.90	1.11	2.10	20.15
	Economic crops	7.55	3.67	2.25	0.44	1.34	92.04
	Other nonvegetated areas	38.44	18.66	4.63	0.22	7.34	21.30
D	Rice	70.31	36.03	0.95	1.05	1.25	13.85
	Trees	32.12	16.46	1.96	0.51	2.20	13.24
	Water bodies	62.61	32.08	1.15	0.87	1.94	7.00
	Economic crops	5.95	3.05	15.12	0.07	4.38	21.55
	Other nonvegetated areas	24.15	12.38	8.49	0.12	18.53	13.50
E	Rice	55.62	31.15	1.11	0.90	1.40	13.64
	Trees	41.63	23.32	2.45	0.41	1.97	14.60
	Water bodies	43.74	24.50	0.78	1.29	2.60	22.36
	Economic crops	5.88	3.29	5.44	0.18	1.58	65.06
	Other nonvegetated areas	31.67	17.74	5.27	0.19	8.29	19.73

3.5. Influence of the Mixed-Pixel

In the five ground-truth sites (Table 5), the average area percentage of water bodies was the highest among the 5 land cover types (larger than 0.71 in five field sites), followed by rice and trees (larger than 0.64 and 0.41, respectively). About 36.03% area of site D was rice, compared with the smallest proportion of 19.66% in site C. The site D had the smallest trees area proportion of 16.46%. The economic crops had the smallest area, and its average patch size in site D was only 0.07 ha, smaller than the area of one pixel of HJ-1 CCD image (30 m × 30 m); the high values of MNN also indicated the highly scattered status of the economic crops (see also Figure 2), and economic crops had the largest value of MNN of the five categories for five field sites. The other nonvegetated areas had the highest AWMSI (except site B), reflecting the complex shape of the road system. The fragmentation statuses of rice, indicated by PD, AWMSI and MNN, of sites D and E were the lowest compared with sites A–C; while the site D and E had the highest MPS. The sites D and E had less fragmented degrees compared with site C, while sites A and B had intermediate level of fragmentation statuses of rice.

Figure 8 showed the ratios of the pixels which were classified as rice field in the HJ CCD images to the total pixel number (ground-truth data, pixels in which the area proportion of rice field is greater than 50%) in each grade for sites A–E. It is obvious that the recognition ratio increased as the area proportion of rice field in pixel ascending, i.e., the grade 75%–100% had the highest classification accuracy. The site D and E had the highest recognition ratio in each grade, while site C the lowest. This result also demonstrated the difficulties in classification in fragmented areas where the mixed-pixel problems were more serious.

More than 69.03% commission pixels concentrated at the boundaries, while at least 63.89% omission pixels lay on the boundaries; that is, most of the misclassified pixels concentrated at the boundaries of the rice fields (Table 6 and Figure 9). The omission pixel numbers of sites D–E were 40 and 36 respectively, obviously less than sites A–C (65, 70 and 80 respectively). The misclassification error was largely determined by the commission error. As shown in Figure 9, the commission error pixels (red color) were more than the omission ones (blue color).

Figure 8. Recognition ratios for the ground-truth sites A–E, compared with the classification results of HJ CCD images.

Figure 8. Recognition ratios for the ground-truth sites A–E, compared with the classification results of HJ CCD images.

Table 6. Statistics of the pixel numbers in classification for sites A–E. The boundary pixels are the pixels which contain certain area of rice field and interface with the other land cover types; the commission pixels represent pixels in which the rice field area proportion is less than <50%, but is misclassified as rice field; and the omission pixels are pixels in which the rice field area proportion is great than 50% but is wrongly classified as the other land cover types.

Statistics of Classification	A	B	C	D	E
Rice pixel number	779	995	678	942	742
Boundary pixel number	476	594	433	542	369
Commission pixel number	150	217	155	144	119
Commission pixels on boundary	122	169	107	101	91
Omission pixel number	65	70	80	44	36
Omission pixels on boundary	47	57	53	40	23
Classification error(%) 1	27.60	28.84	34.66	19.96	20.89

Note: ¹ The classification error is calculated as (omission number + commission number)/rice pixel number.

Table 6. Statistics of the pixel numbers in classification for sites A–E. The boundary pixels are the pixels which contain certain area of rice field and interface with the other land cover types; the commission pixels represent pixels in which the rice field area proportion is less than <50%, but is misclassified as rice field; and the omission pixels are pixels in which the rice field area proportion is great than 50% but is wrongly classified as the other land cover types.

Statistics of Classification	A	B	C	D	E
Rice pixel number	779	995	678	942	742
Boundary pixel number	476	594	433	542	369
Commission pixel number	150	217	155	144	119
Commission pixels on boundary	122	169	107	101	91
Omission pixel number	65	70	80	44	36
Omission pixels on boundary	47	57	53	40	23
Classification error(%) 1	27.60	28.84	34.66	19.96	20.89

Note: ¹ The classification error is calculated as (omission number + commission number)/rice pixel number.

Figure 9. Spatial distribution of the misclassified pixels in sites (A–E).

Figure 9. Spatial distribution of the misclassified pixels in sites (A–E).

4. Discussion

It is generally acknowledged that using a single-temporal image to well discriminate a specific kind of crop at various phenology stages from the other vegetation (or land cover types) is an enormous challenge [19,60,61]. However, using the spectral characteristics (or vegetation indices) determined by the key phenologies of a specific crop species, i.e., multi-temporal remote sensing imageries, is a promising way to improve the classification accuracy [62,63]. To effectively discriminate the rice field in eastern plain region of China, where the rice field is generally fragmented and irregular due to the topography and widely distributed water bodies and road networks, a specifically designed stepwise remote sensing classification strategy was applied in this study.

The time-series EVI2 data for the major land cover types in the study area were built from the HJ-1 A/B CCD imageries and the S-G filters was applied to smooth the EVI2 time-series. With the reference field campaign data, the EVI2 showed efficient discriminating capability in capturing the spectral differences between SCR and the other land cover types during the key SCR phenology stages (Figure 5). It is prominent that the EVI2 of SCR increased rapidly during the transplanting and ear differentiation (including early heading) stages, and the temporal resolution of HJ-1 A/B CCD data was testified to be suitable to capture these features. The stepwise classification algorithm proposed in this study can be seen as an exemplar of the decision tree classification category, and it outperformed the parametric (MLC) and nonparametric (SVM) classification algorithms, respectively (Table 3 and Table 4). By using EVI2 instead of the reflectance data, the classification accuracies improved to certain extents for both of MLC and SVM. The results also implied that by treating the satellite-derived vegetation classification information hierarchically, the mixtures among spectral feature spaces can be effectively alleviated. For MLC and SVM, total six scenes during SCR transplanting to early reproductive stages (from 2012/06/29 to 2012/09/02) were used, including the transplanting, vegetative to reproductive transition phases. However, it is noteworthy that time-series EVI2 of SCR during this period increased rapidly and intersected with the EVI2s of all the other land cover types, except water bodies (Figure 5). Therefore, the classification accuracies of MLC and SVM should unavoidably be decreased, because both of the methods treated the spectral signatures contained in the six scenes collectively.

The influence of the mixed-pixel is a primary concern in remote sensing classification practices. We used five ground-truth sites as an example and analyzed the relationship between the purity of pixels (measured as the area proportion of rice field in a specific cell) and the corresponding recognition ratios. The mixed-pixel analysis showed that the recognition ratio was positively correlated with the rice field area proportion at each ground-truth site (Table 5 and Figure 8). The sites D and E showed the best recognition ratio of SCR among the five ground-truth sites, and it is in accordance with the fragmentation statuses indicated by the landscape indices (Table 5). It is not unexpectedly that as the area proportion of rice field increased in each cell, the possibility of misclassification decreased consequently, especially for the grade of 75%–100% (refer in particular to rice field).

As large part of the classification error can be attributed to the influence of mixed-pixels where the area proportion of rice field was less than 75%, and most of the mixed-pixels concentrated at the boundaries of the rice fields (Table 6 and Figure 9). We further analyzed the classification error caused by the commission and omission errors due to the mixed-pixels and boundary effects, respectively. The results showed that the ratio of the edge pixels to the total rice pixel number correlated with the fragmentation states of each site, i.e., the number of the edge pixel was positively correlated with land fragmentation states of each site due to the increased rice field perimeter. As a consequence, the classification errors of sites D–E were less than sites A–C as shown in Table 6. For rice fields, the misclassification caused by the commission errors was more common, compared with the omission errors (i.e., cells in which rice field area was less than 50% but was classified as rice field, see Figure 8).

However, it should be noted that due to the existence of spatial autocorrelation, the classification accuracy reported in this study may be overestimated [64]. Spatial autocorrelation might be present due to large pixel size [65] or points sampled in close proximity [66]. To avoid the artificially increased classification accuracy caused by the random cross-validation using autocorrelated dataset, more than one permanent training/test dataset should be utilized in accuracy assessment [64]. In this study, five ground-truth sites with different land cover percentages were selected for classification accuracy assessment, however the authors acknowledged that the autocorrelation may still unavoidable and quantitative evaluation of its influence is still a challenge. Further studies should be focused on field data collection, with subsampling and cross-validation like k-fold method [64] to improve the classification accuracy assessment.

The extrapolation of the findings in this study must be cautious due to various changes in the environmental factors (e.g., dry or wet) and vegetation status in different regions and years. The aim of this study was to provide a general methodology in the classification of single-cropped rice. However, when applying it to another region or year, the VI thresholds, which are used to distinguish different land cover types, must be decided according to the specific time series satellite images, i.e., the VI thresholds and the timestamps (according to the key phenologies) are variable.

5. Conclusions

In this study, we applied a simple but robust stepwise algorithm to estimate the single cropped rice (SCR) growing area in irregular and fragmented regions. The multi-temporal HJ-1A/B images and specific signatures of EVI2 at the key phenology stages, i.e., the transplanting and the vegetative to reproductive transition phases, of the SCR were used to classify the rice fields from the other land use types with satisfactory results, compared with the traditional MLC and SVM methods. Due to the fragmented land use composition in the study area, we also assessed the influence of mixed-pixel quantified by using the landscape indices and it showed that the classification accuracy ratio of rice field was positively correlated with its compactness. We showed that by making full use of the key phenological information, and under the support of high-temporal resolution remote sensing data, e.g., HJ-1A/B, the SCR can be mapped at a relative high confidence. The crucial point in the proposed method was the construction of high-quality time-series VI curves, which were then used to identify the key phenology stages to differentiate different land cover types. However, due to the variation of environmental factors and the corresponding changes of vegetation status, due care should be taken when extrapolating the results to other regions or periods. Additionally, we noted that the influence of spatial autocorrelation should also be taken into consideration in classification accuracy evaluation in further study.