Dynamic Mapping of Rice Growth Parameters Using HJ-1 CCD Time Series Data

The high temporal resolution (4-day) charge-coupled device (CCD) cameras onboard small environment and disaster monitoring and forecasting satellites (HJ-1A/B) with 30 m spatial resolution and large swath (700 km) have substantially increased the availability of regional clear sky optical remote sensing data. For the application of dynamic mapping of rice growth parameters, leaf area index (LAI) and aboveground biomass (AGB) were considered as plant growth indicators. The HJ-1 CCD-derived vegetation indices (VIs) showed robust relationships with rice growth parameters. Cumulative VIs showed strong performance for the estimation of total dry AGB. The cross-validation coefficient of determination (RCV) was increased by using two machine learning methods, i.e., a back propagation neural network (BPNN) and a support vector machine (SVM) compared with traditional regression equations of LAI retrieval. The LAI inversion accuracy was further improved by dividing the rice growth period into before and after heading stages. This study demonstrated that continuous rice growth monitoring over time and space at field level can be implemented effectively with HJ-1 CCD 10-day composite data using a combination of proper VIs and regression models.


Introduction
Rice is recognized as one of the most important crops in the world and the staple food for more than half of the world's population [1].Leaf Area Index (LAI) is defined as one-half the total green leaf area per unit of ground surface area [2], and the total weight of plant parts above the ground per unit area at any particular time is called the aboveground biomass (AGB).LAI and AGB are both important vegetation biophysical variables that are widely applied in crop growth monitoring and yield estimation [3].Timely acquisition of crop information on LAI and AGB is vital for the design of dynamic maps of rice growth parameters.Such data could support agronomic decisions that ensure plant health, yield stability, and optimization of economic benefits.
Traditional crop growth monitoring methods are limited in space and time, especially in meeting the requirements of continuous crop condition monitoring at regional scales.Remote sensing techniques offer a feasible tool for parameter regionalization.Lu made a review of the potential and challenge of remote sensing-based forest biomass estimation; Koppe et al. did a study on the estimation of winter wheat biomass using 30 m resolution data from Hyperion on regional scale; and Liu et al. explored the uncertainties in green LAI estimation from Landsat images over multiple growing seasons with several crop types [4][5][6].However, few studies have been conducted on the dynamic mapping of crop growth parameters, which is important for practical applications.Currently, the most widely used method for estimating LAI or AGB of crops from remote sensing data is through the use of statistical relationships between the parameters and spectral vegetation indices (VIs) derived from many different kinds of optical remote sensing data [7][8][9] or other nonparametric methods [10,11].Some researchers have explored the physical model inversion methods with satisfactory results [10,12].These inversion models usually contain many parameters and often require extensive details of meteorological and/or agricultural knowledge to run, making real-time inversion on a regional scale unsuitable.However, less is known about the application of machine learning methods and their potentials in inverting crop biophysical parameters using remote sensing data.
The Advanced Very High Resolution Radiometer (AVHRR), Moderate Resolution Imaging Spectroradiometer (MODIS), Satellite pour l'Observation de la Terre VEGETATION (SPOT-VGT) and Landsat-MSS have demonstrated advantages in rice monitoring from regional to global scales due to wider coverage and relatively longer data archive [13][14][15].Some coarse resolution data such as AVHRR-NDVI, MODIS-NDVI, and MODIS-EVI have been preprocessed and provide standard VI products at a global scale to end users [16][17][18].However, coarse resolution satellite data are generally prone to classification uncertainties as mixed land cover types are difficult to delineate, especially in the highly fragmented rice fields in the eastern plains of China [19].Moreover, in practice, crop LAI was found to be significantly underestimated by the MODIS LAI product [20,21], which makes the implementation of suitable field management practices at a local scale more complicated.
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 [22].However, their cost and relatively long revisit cycles partially offset their advantages in spatial resolution.Specifically, most of the rice fields located in the subtropical regions of China are influenced by the monsoon climate, leading to inadequate acquisitions of cloud-free remote sensing data during monsoon weather conditions.In precision agriculture, where rice growth information is critically needed, high-resolution multi-temporal remote sensing data are advantageous, as the spatial variability captured by remote sensing data is useful for adjusting crop properties while taking into account the local conditions [5].
Two small environment and disaster monitoring and forecasting satellites (HJ-1A/B) were launched by the China Center for Resources Satellite Data and Applications (CRESDA) in 2008.The charge-coupled device (CCD) cameras of these satellites have a swath width of 700 km, four spectral bands with a spatial resolution of 30 m in the visible bands, and a revisit cycle of four days (the revisit cycle of the constellation is two days).This fine spatial resolution and short repeat cycle facilitate the availability of images at all critical growth stages of crops [23,24].
Vegetation Indices (VIs) are robust empirical measures of vegetation activity at the land surface.They are designed to enhance the vegetation signal from measured spectral responses by combining two (or more) different wavebands, often in the red (0.6-0.7 µm) and NIR wavelengths (0.7-1.1 µm), and are widely used in the studies of crop land classification, yield prediction, phenology, and crop growth monitoring [25][26][27].The normalized difference vegetation index (NDVI) [28] is by far the most widely used index in the literature, and is advantageous for studying historical changes; however, it is sensitive to canopy background variations and saturates in relatively high-vegetated areas [29,30].The enhanced vegetation index (EVI) [30] and the two-band enhanced vegetation index (EVI2) [31] differ from NDVI by attempting to correct for atmospheric and background perturbation, and appear to be superior in discriminating subtle differences in areas of high vegetation density [32].Cumulative VIs have been proven as a surrogate for absorbed photosynthetically active radiation (APAR), which is proportional to total biomass [33,34].
The objectives of this study were to evaluate the use of vegetation indices for continuous rice growth monitoring over multiple years and at a regional scale using a newly constructed HJ-1 CCD 10-day composite data, explore the traditional statistical models and nonparametric methods for LAI and AGB retrieval based on extensive field campaign datasets, and propose a dynamic mapping method for rice growth for in situ application in agricultural practice.

Study Site
The study was conducted in Deqing County, which lies in the western part of the Hangjiahu Plain, Zhejiang Province, Southeast China (Figure 1).This area is characterized by a tropical monsoon climate.The mean annual temperature ranges between 13 • C and 16 • C and there is an annual total precipitation of 1379 mm.Deqing has a total area of 936 km 2 , and according to statistical data from the local agriculture department, more than 91% of the crop area is single-cropped rice (SCR).However, due to the countless lakes, ponds, and rivers scattered throughout this region, and the well-developed road networks, an irregular crop land phenomenon is typical in the study area [19].
Remote Sens. 2016, 8, 931 3 of 18 10-day composite data, explore the traditional statistical models and nonparametric methods for LAI and AGB retrieval based on extensive field campaign datasets, and propose a dynamic mapping method for rice growth for in situ application in agricultural practice.

Study Site
The study was conducted in Deqing County, which lies in the western part of the Hangjiahu Plain, Zhejiang Province, Southeast China (Figure 1).This area is characterized by a tropical monsoon climate.The mean annual temperature ranges between 13 °C and 16 °C and there is an annual total precipitation of 1379 mm.Deqing has a total area of 936 km 2 , and according to statistical data from the local agriculture department, more than 91% of the crop area is single-cropped rice (SCR).However, due to the countless lakes, ponds, and rivers scattered throughout this region, and the well-developed road networks, an irregular crop land phenomenon is typical in the study area [19].

Measurement of Crop Parameters
Single-cropped rice is the most abundant crop in Deqing County.Most small patches of rice fields were cultivated by small-scale individual farmers, whereas relatively larger fields were managed by large-scale rice farmers, leading to time inconsistency in rice cultivation and management measures, including different transplanting dates and different rice varieties.Extensive crop parameters experiments were conducted during the rice growing seasons (June to November) for the years 2012 and 2013.Field and laboratory measurements included LAI, AGB, and plant density.The corresponding rice phenology dates and sample numbers for each date were recorded as well (Table 1).
Rice is cultivated on flooded soils.At the early stages of growth, rice fields are a mixture of green rice plants and open water [35].Water depth generally varies from 2 to 15 cm.As rice continues to grow, new leaves emerge.The development of tillers is slow at first and becomes faster until the maximum tiller number is attained at the highest plant density.About 50 to 60 days after transplanting or sowing (the transplanting or sowing dates for 2012 and 2013 were 21 June and 18 June, respectively), rice plant canopies cover most of the soil area.The leaves continue to increase in total area and generally remain green until the heading stage.Due to intraspecific competition at the reproductive stage, the tiller number drops and the rice leaves gradually become yellow.This stage is also associated with an increase in the dry aboveground biomass of rice until maturity (Figure 2) [36].

Measurement of Crop Parameters
Single-cropped rice is the most abundant crop in Deqing County.Most small patches of rice fields were cultivated by small-scale individual farmers, whereas relatively larger fields were managed by large-scale rice farmers, leading to time inconsistency in rice cultivation and management measures, including different transplanting dates and different rice varieties.Extensive crop parameters experiments were conducted during the rice growing seasons (June to November) for the years 2012 and 2013.Field and laboratory measurements included LAI, AGB, and plant density.The corresponding rice phenology dates and sample numbers for each date were recorded as well (Table 1).
Rice is cultivated on flooded soils.At the early stages of growth, rice fields are a mixture of green rice plants and open water [35].Water depth generally varies from 2 to 15 cm.As rice continues to grow, new leaves emerge.The development of tillers is slow at first and becomes faster until the maximum tiller number is attained at the highest plant density.About 50 to 60 days after transplanting or sowing (the transplanting or sowing dates for 2012 and 2013 were 21 June and 18 June, respectively), rice plant canopies cover most of the soil area.The leaves continue to increase in total area and generally remain green until the heading stage.Due to intraspecific competition at the reproductive stage, the tiller number drops and the rice leaves gradually become yellow.This stage is also associated with an increase in the dry aboveground biomass of rice until maturity (Figure 2) [36].Fourteen rice sample plots with areas larger than 200 × 200 m 2 in Deqing County were randomly chosen in order to measure the vegetation characteristics (Figure 1).A handheld global positioning system (GPS) receiver (Trimble Juno-SB, Trimble Navigation Ltd., Sunnyvale, CA, USA) was used to record the precise location and geometry of every sample plot with a positioning accuracy of approximately ±5 m.At every measurement time, five quadrants of 0.5 m × 0.5 m were selected from each plot, and all single-cropped plants within selected quadrants were harvested (Figure 3).Fourteen rice sample plots with areas larger than 200 × 200 m 2 in Deqing County were randomly chosen in order to measure the vegetation characteristics (Figure 1).A handheld global positioning system (GPS) receiver (Trimble Juno-SB, Trimble Navigation Ltd., Sunnyvale, CA, USA) was used to record the precise location and geometry of every sample plot with a positioning accuracy of approximately ±5 m.At every measurement time, five quadrants of 0.5 m × 0.5 m were selected from each plot, and all single-cropped plants within selected quadrants were harvested (Figure 3).* The heading stage of single-cropped rice (SCR) according to the field campaign.LAI: leaf area index; AGB: aboveground biomass.The aboveground parts of rice plants were separated into leaves, stems, and ear (after ear differentiation) in the lab.The leaf areas of single-cropped rice were measured using a leaf area meter (LI-3000C, with conveyor belt assembly, LI-3050C; Li-Cor, Inc., Lincoln, NE, USA).All plant parts were weighted by an electronic scale in fresh conditions at first.Then they were put into the oven at 105 °C for 30 min and kept at 70 °C for 72 h until they reached constant weight.The dry aboveground biomass was weighted with the same electronic balance and scaled to total dry AGB m 2 using plant density, which was measured as the average number of plants of five quadrats in each sample plot.The average LAI and AGB of five quadrats in each rice sample plot represent the vegetation characteristic values of that sample plot.

Remote Sensing Data and Vegetation Indices (VIs)
HJ-1 CCD data from June to November in 2012 and 2013 over the study area were downloaded from the China Centre for Resources Satellite Data and Application (CRESDA).The sensor The aboveground parts of rice plants were separated into leaves, stems, and ear (after ear differentiation) in the lab.The leaf areas of single-cropped rice were measured using a leaf area meter (LI-3000C, with conveyor belt assembly, LI-3050C; Li-Cor, Inc., Lincoln, NE, USA).All plant parts were weighted by an electronic scale in fresh conditions at first.Then they were put into the oven at 105 • C for 30 min and kept at 70 • C for 72 h until they reached constant weight.The dry aboveground biomass was weighted with the same electronic balance and scaled to total dry AGB m 2 using plant density, which was measured as the average number of plants of five quadrats in each sample plot.The average LAI and AGB of five quadrats in each rice sample plot represent the vegetation characteristic values of that sample plot.

Remote Sensing Data and Vegetation Indices (VIs)
HJ-1 CCD data from June to November in 2012 and 2013 over the study area were downloaded from the China Centre for Resources Satellite Data and Application (CRESDA).The sensor characteristics are presented in Table 2.A total of 100 scenes of HJ-1A/B images (46 scenes and 54 scenes in 2012 and 2013, respectively) with cloud cover less than 20% during the single-cropped rice phenology in the study area were selected for the following procedures.All the images collected were processed through radiometric calibration, atmospheric correction, and geometric correction.The formulas and coefficients for radiometric calibration were collected from the raw image package.The atmospheric correction was performed using the Moderate Resolution Transmission (MODTRAN) 4 model integrated in the Fast Line-of-sight Atmospheric Analysis of Spectral Hypercubes (FLAASH) module in the Environment for Visualizing Images (ENVI) package.The images were co-registered using the Second National Soil Survey Vector Map (scale 1:10,000) provided by the Deqing County Land and Resources Bureau, and the geometric accuracy was less than 0.5 pixel (15 m).Only cloud-free images that corresponded to each field campaign date within three days were selected for LAI modeling.Using Savitzky-Golay (S-G) filters, daily VI time series were created from all the HJ-1 CCD images (100 scenes in two years).These images were used to verify the AGB estimation before dynamic mapping [23,37].After modeling and model validation, all of the images were processed with cloud detection and maximum-value composite (MVC) methods to produce a new HJ-1 CCD 10-day composite data for use in the timely dynamic mapping of rice growth.Cloud-free pixels usually have lower reflectances between 0 and 0.3 in the red band, whereas thick cloud-contaminated pixels have higher reflectances distributed over large ranges between 0 and 1.7 [38].The threshold value for HJ-1 CCD data was selected as 0.25 for the red band after statistics calculated from 1214 cloud pixels and 1029 cloud-free pixels randomly selected from HJ-1 CCD images.The flagged cloudy pixels were set to no data.MVC was initially developed for AVHRR and MODIS VIs data and has been widely used to minimize atmospheric effects, including residual clouds [39,40].MVC was designed based on the assumption that observations with a higher component of haze or cloud will show a higher reflectance in the red band, thus a lower VI value to distinguish the edge between land and residual cloud [38,41].For modeling of LAI from HJ-1 CCD, there were 191 matching LAI sample sites and 194 matching AGB sample sites.
In addition to the selected VIs, NDVI, and EVI2, cumulative VIs were adopted for the estimation of AGB.
The matched cumulative VIs were calculated from day of year (DOY) 173 and 169 for 2012 and 2013, respectively, approximately representing the transplanting date or sowing date of SCR.

Deriving LAI and AGB via VIs
Pixel-based VI values from the corresponding HJ-1 CCD images at the 10 field measurement plots were extracted from the images.As the SCR LAI drops after heading, the growth stages of SCR were divided into before heading and after heading in order to improve the estimation results.
Best fit linear and non-linear relationships among VIs, LAI, and AGB were established first.In addition to the five traditional regression equations (linear, exponential, power, logarithmic, and quadratic polynomial regression), two machine learning methods, i.e., a back propagation neural network (BPNN) and a support vector machine (SVM), were applied in model construction.The neural network approach has the advantages of nonlinearity, input-output mapping, adaptivity, generalization, and fault tolerance [27].A SVM for regression analysis is accomplished by solving a convex optimization problem, more specifically a quadratic programming problem [42].BPNN and SVM have demonstrated their abilities in dealing with complicated datasets in many fields like classification [43,44] and data mining [45,46].
In this work, the neural network function in the MATLAB package was used for the BPNN regression; LIBSVM 3.20 [47] and LIBSVM-FarutoUltimate toolboxes (software available at http://www.ilovematlab.cn)were used for the SVM analyses.The kernel function used in SVM was the radial basis function (RBF) kernel [42].The two tuning parameters used in the RBF kernel for LIBSVM are "cost" (C) and "gamma" (γ).The ranges of γ and C were set to [-15, 15] and [-15, 15] for regression, and the steps were 0.5 and 0.5, respectively.To find the optimal parameter, an iterative grid search was applied in which the evaluation was provided by a 10-fold cross-validation for both BPNN and SVM regression.
The leave-one-out cross-validation (LOO-CV) method was implemented to test the prediction capability of the models.This method sets a single observation from the original sample as validation data, and the remaining observations as training data.The performances of the models mentioned above were compared on the basis of cross-validation coefficient of determination (R 2 CV ) and relative cross-validation root mean square error (RRMSE CV ).RRMSE CV facilitates comparison of the accuracy of each estimation model and was calculated using the following equation: where n is the number of samples and P i and Q i are the predicted and observed values, respectively.Q is the observed mean value.All analyses were conducted in the MATLAB program environment.Furthermore, we verified the relationship between the 10-day composite VI data with AGB.The AGB curve can be described with logistic curves [36]; the particular form of the logistic equation available was written as: where W is the dry aboveground biomass weight of the crop, t is time, k is the growth rate, A + C is the maximum theoretical AGB weight, and b is a constant term.The accuracy of the model result was determined by R 2 and was accomplished using the MATLAB program.The cumulative VIs were calculated from all 10-day composite data and the corresponding AGB values were estimated from the aforementioned logistic function.The relationships between daily cumulative VIs and 10-day composite cumulative VIs were analyzed, and the best regression equation was built for the dynamic mapping application.

Dynamic Mapping
The spatial distribution and temporal dynamics of LAI and AGB of SCR were determined using the best regression model based on the HJ-1 CCD 10-day composite data during the SCR growth season.Before the dynamic mapping process, the planting area of single-cropped rice was estimated using a stepwise classification strategy proposed by Wang et al. [19].The total acreage of the single-cropped rice was about 94.0 km 2 .The levels of SCR growth parameters were visualized with different colors.

Relationships between VIs and Rice Growth Parameters
To quantify the effects of NDVI, EVI2, cumulative NDVI, and cumulative EVI2 on SCR LAI and AGB, we extracted the value of remote sensing data corresponding to the field campaign date.The correlation between rice growth parameters and each feature before and after the heading stage are analyzed in Table 3.The VIs were significantly correlated with LAI and AGB, except that between VIs and AGB at the all-growth stage.Dividing the SCR growth period into before heading and after heading, the correlation coefficients between LAI and image features increased to a great extent.EVI2 showed the highest correlation coefficients with LAI at the before heading stage (i.e., 0.86).NDVI had the best relationship with LAI at the after heading stage and its correlation coefficient was 0.67, which was higher than EVI2.NDVI and EVI2 did not show a stable, positive correlation with AGB.Among the first two spectral features in Table 3, EVI2 showed the best correlation coefficient of 0.83 with AGB at the before heading stage.Even though cumulative NDVI and cumulative EVI2 had an overall high correlation at all stages, their correlation coefficient, taking the after heading stage alone, was relatively low.

LAI and AGB Regression Model Analysis
LAI and AGB regression models were individually constructed using five traditional regression equations (linear, exponential, power, logarithmic, and quadratic polynomial regression) and two machine learning methods (BPNN and SVM) with VIs.Cumulative EVI2 and cumulative NDVI were only used for AGB regression.The results of traditional models and two machine learning methods are listed in Table 4.
For the LAI regression analysis during all the growth stages of SCR, traditional statistical models did not perform as well as BPNN and SVM.The SVM model constructed the best regression function in describing this group; the R 2 CV was 0.47 and RRMSE CV was 10.19 when using NDVI as the input dataset.The cross-validation coefficient of determination became higher after dividing the growth period into before heading and after heading groups based on the field campaign data.By dividing the heading data, the growth phases of SCR can be classified into vegetative growth and reproductive growth.EVI2 were superior to NDVI among all the regression methods during the before heading stage.The BPNN model was recommended in LAI estimation.After heading, rice leaves begin to fall due to etiolation and senescence.The SVM-based regression function could explain most of the variance of the dependent variable and the RRMSE CV was also the lowest, i.e., 7.08.EVI2, NDVI, cumulative EVI2, and cumulative NDVI were used in the SCR dry AGB estimation.However, the correlation of EVI2 and NDVI with AGB was not significant across the growth stage.Cumulative VIs had a satisfactory result with AGB, especially during all stages and the before heading stage.The best traditional regression method was the quadratic polynomial, with R 2 CV values of 0.93 at all stages and 0.92 at the before heading stage, using cumulative NDVI as the independent variable.BPNN and SVM showed great abilities in AGB regression as well.SVM could get the best result at the after heading stage.
Based on the above analysis, the best strategy for estimating LAI is to divide the SCR growth period into before heading and after heading stages.EVI2-BPNN regression was selected for the before heading LAI estimation, and NDVI-SVM was selected for the after heading LAI estimation.The cumulative NDVI-based quadratic polynomial fit function was adopted for the prediction of AGB at all stages (Figure 4a-c).These changes have no material impact on the conclusions of our paper.We apologize to our readers for the inconvenience.The manuscript will be updated and the original will remain online on the article webpage.These changes have no material impact on the conclusions of our paper.We apologize to our readers for the inconvenience.The manuscript will be updated and the original will remain online on the article webpage.

Dynamic Mapping Method of Rice Growth Monitoring
Considering the dynamic mapping approach, the HJ-1 CCD images were cloud detected and a per-pixel VI maximum-value composite generated to form a new 10-day composite data (three data points per month) from June to November 2012 and 2013.As the daily cumulative NDVI is difficult to obtained, the accuracy of 10-day composite data based on cumulative NDVI is crucial.All the time series measured AGB values per sample were fitted with logistic curves.The results are shown in Figure 5.
According to the time series AGB values per sample and the 10-day composite data derived from cumulative NDVI, the regression equation was acquired using the quadratic polynomial fit function; the specific regression models are presented in Figure 6.

Dynamic Mapping Method of Rice Growth Monitoring
Considering the dynamic mapping approach, the HJ-1 CCD images were cloud detected and a per-pixel VI maximum-value composite generated to form a new 10-day composite data (three data points per month) from June to November 2012 and 2013.As the daily cumulative NDVI is difficult to obtained, the accuracy of 10-day composite data based on cumulative NDVI is crucial.All the time series measured AGB values per sample were fitted with logistic curves.The results are shown in Figure 5.
According to the time series AGB values per sample and the 10-day composite data derived from cumulative NDVI, the regression equation was acquired using the quadratic polynomial fit function; the specific regression models are presented in Figure 6.The dynamic maps of SCR growth conditions in Deqing County are presented in Figures 7-10.The composite VI images represent the highest value of each pixel during the 10 days, and estimated the SCR growth conditions three times a month for field level practices.The dynamic mapping can realize the near real-time observation at a regional scale with satisfactory precision.At different years, there is an observed variation in modeled LAI values at the corresponding time.For example,

Dynamic Mapping Method of Rice Growth Monitoring
Considering the dynamic mapping approach, the HJ-1 CCD images were cloud detected and a per-pixel VI maximum-value composite generated to form a new 10-day composite data (three data points per month) from June to November 2012 and 2013.As the daily cumulative NDVI is difficult to obtained, the accuracy of 10-day composite data based on cumulative NDVI is crucial.All the time series measured AGB values per sample were fitted with logistic curves.The results are shown in Figure 5.
According to the time series AGB values per sample and the 10-day composite data derived from cumulative NDVI, the regression equation was acquired using the quadratic polynomial fit function; the specific regression models are presented in Figure 6.The dynamic maps of SCR growth conditions in Deqing County are presented in Figures 7-10.The composite VI images represent the highest value of each pixel during the 10 days, and estimated the SCR growth conditions three times a month for field level practices.The dynamic mapping can realize the near real-time observation at a regional scale with satisfactory precision.At different years, there is an observed variation in modeled LAI values at the corresponding time.For example, The dynamic maps of SCR growth conditions in Deqing County are presented in Figures 7-10.The composite VI images represent the highest value of each pixel during the 10 days, and estimated the SCR growth conditions three times a month for field level practices.The dynamic mapping can realize the near real-time observation at a regional scale with satisfactory precision.At different years, there is an observed variation in modeled LAI values at the corresponding time.For example, in mid-July 2013, SCR showed a higher LAI value compared with that observed at the same time in 2012 (Figures 7b and  8b), which may be caused by slight differences in the transplanting or sowing dates (compared with 2012).Both years showed a peak LAI value at the start of the reproductive stage, and after the heading stage the LAI dropped due to etiolation and leaf senescence.In mid-November, many fields were harvested and LAI values were relatively low (Figure 7m,n and Figure 8m,n).The dry AGB of rice showed continuous increase during the growth period (Figures 9 and 10).Different colors represent different ranges of AGB values, and could provide information for field managers and policymakers at regional scales.The overall system of dynamic mapping is shown in Figure 11.The specific processes were integrated in the MATLAB program.
Remote Sens. 2016, 8, 931 11 of 18 in mid-July 2013, SCR showed a higher LAI value compared with that observed at the same time in 2012 (Figures 7b and 8b), which may be caused by slight differences in the transplanting or sowing dates (compared with 2012).Both years showed a peak LAI value at the start of the reproductive stage, and after the heading stage the LAI dropped due to etiolation and leaf senescence.In mid-November, many fields were harvested and LAI values were relatively low (Figures 7m,n and 8m,n).The dry AGB of rice showed continuous increase during the growth period (Figures 9  and 10).Different colors represent different ranges of AGB values, and could provide information for field managers and policymakers at regional scales.The overall system of dynamic mapping is shown in Figure 11.The specific processes were integrated in the MATLAB program.

Discussion
The LAI remote sensing products now available are usually of coarse spatial resolution and LAI values have mostly been underestimated due to the mixed-pixel problem [48].Little research has been done on the remote sensing of rice AGB at higher spatio-temporal scales and multi-season dynamic mapping of the entire growth process.Considering the growing importance of biophysical parameters in shaping local agronomic activities, there is an increasing need for the precise monitoring of rice crop growth conditions, especially in the eastern plains region of China, where rice fields are relatively irregular and fragmented.The current study is focused on the near-real-time dynamic mapping of rice growth using multi-temporal middle-to high-resolution remote sensing time series data for regional-scale applications.A comprehensive field campaign was carried out to record the growth parameters of single-cropped rice, including LAI, AGB, and plant density, and was used in model construction and verification.VIs were derived from HJ-1 CCD images for the purpose of growth parameters inversion.
Time series NDVI and EVI2 were calculated from HJ-1 CCD images, which corresponded to the field campaigns within three days.Vegetation indices combined with reflectance bands showed a high correlation with LAI and AGB (Table 3).VIs have the advantage of obtaining relevant information rapidly and easily, and perform better than a single spectral band in crop monitoring [24,49,50].In situ measured LAI reached a maximum, while in situ measured dry AGB continued to increase, and the AGB includes more photosynthetically inactive components (e.g., stem, ear) than the leaf area, which is likely to affect the relationship between VIs and the photosynthetically active components [34].In this situation, cumulative VIs were proposed as a better approach for the estimation of AGB.
The regression model chosen for the derivation of SCR growth parameters is another important issue in model construction.The traditional models (linear, exponential, power, logarithm, and quadratic polynomial regression) are commonly used in crop LAI and dry AGB estimation in previous studies [25,51].Partial least squares regression and stepwise multiple regression models have been tested as feasible tools in information extraction [52], but at the same time introduce a lot of features into the model that result in difficulties in application.There is increasing interest in data mining technology in remote sensing, which is considered to be a useful tool to unravel the non-linear relationship between canopy spectral reflectance and crop growth conditions [11,42].However, few studies have been performed on the application of the technology in agriculture.Backpropagation neural networks (BPNN) and support vector machines (SVM) were applied in model construction in our research.We found that BPNN and SVM performed better than the traditional models in LAI estimation, demonstrating that machine learning methods are potentially useful to understand and predict optical interactions over a wide range of rice canopy LAI.Although the quadratic polynomial regression showed good performance in AGB estimation, the BPNN and SVM methods can also reach comparably high accuracy.This demonstrates that the non-linear methods are efficient in establishing relationships between remote sensing data and rice biophysical parameters.

Discussion
The LAI remote sensing products now available are usually of coarse spatial resolution and LAI values have mostly been underestimated due to the mixed-pixel problem [48].Little research has been done on the remote sensing of rice AGB at higher spatio-temporal scales and multi-season dynamic mapping of the entire growth process.Considering the growing importance of biophysical parameters in shaping local agronomic activities, there is an increasing need for the precise monitoring of rice crop growth conditions, especially in the eastern plains region of China, where rice fields are relatively irregular and fragmented.The current study is focused on the near-real-time dynamic mapping of rice growth using multi-temporal middle-to high-resolution remote sensing time series data for regional-scale applications.A comprehensive field campaign was carried out to record the growth parameters of single-cropped rice, including LAI, AGB, and plant density, and was used in model construction and verification.VIs were derived from HJ-1 CCD images for the purpose of growth parameters inversion.
Time series NDVI and EVI2 were calculated from HJ-1 CCD images, which corresponded to the field campaigns within three days.Vegetation indices combined with reflectance bands showed a high correlation with LAI and AGB (Table 3).VIs have the advantage of obtaining relevant information rapidly and easily, and perform better than a single spectral band in crop monitoring [24,49,50].In situ measured LAI reached a maximum, while in situ measured dry AGB continued to increase, and the AGB includes more photosynthetically inactive components (e.g., stem, ear) than the leaf area, which is likely to affect the relationship between VIs and the photosynthetically active components [34].In this situation, cumulative VIs were proposed as a better approach for the estimation of AGB.
The regression model chosen for the derivation of SCR growth parameters is another important issue in model construction.The traditional models (linear, exponential, power, logarithm, and quadratic polynomial regression) are commonly used in crop LAI and dry AGB estimation in previous studies [25,51].Partial least squares regression and stepwise multiple regression models have been tested as feasible tools in information extraction [52], but at the same time introduce a lot of features into the model that result in difficulties in application.There is increasing interest in data mining technology in remote sensing, which is considered to be a useful tool to unravel the non-linear relationship between canopy spectral reflectance and crop growth conditions [11,42].However, few studies have been performed on the application of the technology in agriculture.Backpropagation neural networks (BPNN) and support vector machines (SVM) were applied in model construction in our research.We found that BPNN and SVM performed better than the traditional models in LAI estimation, demonstrating that machine learning methods are potentially useful to understand and predict optical interactions over a wide range of rice canopy LAI.Although the quadratic polynomial regression showed good performance in AGB estimation, the BPNN and SVM methods can also reach comparably high accuracy.This demonstrates that the non-linear methods are efficient in establishing relationships between remote sensing data and rice biophysical parameters.
Most previous studies laid emphasis on rice growth parameters at a single stage, lacking the exploration of consistent growth stages [53,54], and, hence, growth monitoring systems were difficult to find in practice.In order to improve the regression precision, we made use of the specific growth curve of SCR.After the heading stage, as LAI drops due to etiolation and senescence of the leaves, the whole growth period was divided into before heading and after heading stages.The accuracy of LAI estimation was improved with this strategy.EVI2 had better results at the before heading stage than NDVI.This may be caused by the saturation of NDVI in relatively low-vegetated areas; this result agrees with that of Shang et al. [51], who also showed that EVI2 was more resistant to saturation at high biomass range than NDVI.
The cloudy pixels were flagged by a threshold value of 0.25 for the red band of HJ-1 CCD data.The MVC method can further help in the removal of contaminated pixels [55,56].New 10-day composite HJ-1 CCD data were then acquired three times a month for estimating local SCR growth parameters.In the rice growth monitoring and dynamic mapping system, EVI2-BPNN regression is selected for before heading LAI estimation; NDVI-SVM is selected for after heading LAI estimation; cumulative NDVI based on the quadratic polynomial fit function was adopted for the all-stage AGB prediction.The R 2 CV values were 0.93, 0.66, and 0.93, respectively.With this method, rice growth conditions can be visually captured at the regional scale, which will help users compare the overall condition of crops through different growth stages within one year or different years in the same growth stage.
However, it should be noted that the empirical equations built between remote sensing derived VIs and rice growth parameters were partially determined by field experiment data and remote sensing image preprocessing.As the rice planting area may change annually, precise extraction of the crop planting area is crucial for the mapping of growth parameters.There are, therefore, inevitable uncertainties in estimation accuracy.When extrapolating this dynamic mapping system to other regions, considering the variations in field experiment methods and remote sensing data characteristics, necessary verification steps must be taken.According to Koppe et al. and Jin et al., TerraSAR-X (Germany) data and RADARSAR-2 (Canada) data have been successfully used in crop monitoring [54,57].Synthetic aperture radars (SARs) have some advantages for monitoring crop growth status as they are not influenced by the presence of clouds or haze.However, SAR images are usually limited by image processing techniques [25].The method proposed in this research made full use of HJ-1 CCD images with high temporal resolution, acquired concurrently with the field data.Considering that HJ-1 data mostly cover China and some surrounding Asian countries and not the whole globe, the sensors only have four wavebands (Table 2); future studies should be directed towards exploring the potential of SAR (Sentinel-1, ENVISAT-ASAR, RADARSAT-1, etc.) and optical remote sensing data (Gaofen series satellite, FORMOSAT-2, Landsat 8, etc.) for monitoring rice growth parameters other than LAI and AGB.

Conclusions
In this study, we proposed a dynamic mapping method for single-cropped rice growth parameters at regional, annual, and inter-annual scales.The method was developed based on field campaigns conducted in two consecutive rice growing seasons (2012 and 2013), each spanning June to November in Deqing County, Southeast China.Ten-day HJ-1 CCD image composites, giving a total of three image composites per month, were generated for growth parameters retrieval.The field measurements included leaf area index (LAI), aboveground biomass (AGB), and plant density, and were strictly controlled during measurement.The LAI empirical equations were established by dividing the whole SCR growth period into before heading and after heading stages.EVI2 performed better at the fast growth stage (before heading) whereas NDVI showed better accuracy at the after heading stage (relatively slower growth).Two machine learning methods, i.e., backpropagation neural network (BPNN) and support vector machine (SVM), were applied in LAI model construction.Cumulative VIs were proven suitable for the estimation of AGB.Cumulative NDVI based on the quadratic polynomial fit function was adopted for the all-stage AGB prediction.
This study demonstrated the potential of using HJ-1 CCD images in rice growth monitoring.Machine learning methods provided a useful exploratory tool for improvement in the relationships between different combinations of reflectance and crop variables.Much work remains to be done to scale these growth condition estimation relationships across a variety of canopies, and more sources of remote sensing data should be included in the system.Only then will the approach gain sufficient robustness and reliability to be applied at larger scales of agricultural remote sensing.

Figure 1 .
Figure 1.HJ-1A charge-coupled device (CCD) false-color composite image of the study area of Deqing County, Zhejiang Province, East China.

Figure 1 .
Figure 1.HJ-1A charge-coupled device (CCD) false-color composite image of the study area of Deqing County, Zhejiang Province, East China.

Figure 2 .
Figure 2. The single-cropped rice calendar in Deqing County, Zhejiang Province, China: (a) machinery transplants rice seedlings, 21 June 2012; (b) tillering period, about two weeks after transplanting; (c) the tillering period reached the maximum tiller number on 30 July 2012; (d) ear differentiation period, 15 August 2012; (e) heading period, flowering, 31 August 2012; (f) ripening stage, ready to be harvested, 18 November 2012.Photos (a-f) were taken at the same rice sample plot, as there were several rice varieties planted, but the start date of each period is different.

Figure 2 .
Figure 2. The single-cropped rice calendar in Deqing County, Zhejiang Province, China: (a) machinery transplants rice seedlings, 21 June 2012; (b) tillering period, about two weeks after transplanting; (c) the tillering period reached the maximum tiller number on 30 July 2012; (d) ear differentiation period, 15 August 2012; (e) heading period, flowering, 31 August 2012; (f) ripening stage, ready to be harvested, 18 November 2012.Photos (a-f) were taken at the same rice sample plot, as there were several rice varieties planted, but the start date of each period is different.

Figure 3 .
Figure 3. Location of sample quadrats in a sample plot.

Figure 3 .
Figure 3. Location of sample quadrats in a sample plot.

Figure 4 .
Figure 4. Relationships between measured rice leaf area index (m 2 /m 2 ) and dry aboveground biomass (g/m 2 ) at different rice growth stages with VIs.(a) Before heading LAI estimation using EVI2-BPNN regression; (b) after heading LAI estimation using NDVI-SVM regression; (c) all-growth stage AGB estimation using daily cumulative NDVI and based on the quadratic polynomial fit function; (d) all-growth stage AGB estimation using 10-day composite data and based on the cumulative NDVI quadratic polynomial fit function.The black dash lines are the 45° lines, and the red solid lines are the linear regression trend lines.

Figure 4 .
Figure 4. Relationships between measured rice leaf area index (m 2 /m 2 ) and dry aboveground biomass (g/m 2 ) at different rice growth stages with VIs.(a) Before heading LAI estimation using EVI2-BPNN regression; (b) after heading LAI estimation using NDVI-SVM regression; (c) all-growth stage AGB estimation using daily cumulative NDVI and based on the quadratic polynomial fit function; (d) all-growth stage AGB estimation using 10-day composite data and based on the cumulative NDVI quadratic polynomial fit function.The black dash lines are the 45 • lines, and the red solid lines are the linear regression trend lines.

Figure 4 .
Figure 4. Relationships between measured rice leaf area index (m 2 /m 2 ) and dry aboveground biomass (g/m 2 ) at different rice growth stages with VIs.(a) Before heading LAI estimation using EVI2-BPNN regression; (b) after heading LAI estimation using NDVI-SVM regression; (c) all-growth stage AGB estimation using daily cumulative NDVI and based on the quadratic polynomial fit function; (d) all-growth stage AGB estimation using 10-day composite data and based on the cumulative NDVI quadratic polynomial fit function.The black dash lines are the 45 • lines, and the red solid lines are the linear regression trend lines.

Figure 5 .
Figure 5. Dynamic curves of rice dry aboveground biomass (g/m 2 ) based on logarithmic equations: (a) logarithmic curve of one sample plot in 2013; (b) comparison of measured aboveground biomass with estimated values of all the sample plots during 2012-2013 (n = 194).The black curve is the best-fit function for in situ AGB estimation, the black dashed line is the 45° line, and the red solid line is the linear regression trend line.

Figure 6 .
Figure 6.Regression analysis of HJ-1 CCD cumulative VIs versus dry aboveground biomass (g/m 2 ) using a quadratic polynomial fit function: (a) daily cumulative NDVI function and (b) 10-day composite data based on cumulative NDVI function.The blue solid lines are the regression trend lines.

Figure 5 .
Figure 5. Dynamic curves of rice dry aboveground biomass (g/m 2 ) based on logarithmic equations: (a) logarithmic curve of one sample plot in 2013; (b) comparison of measured aboveground biomass with estimated values of all the sample plots during 2012-2013 (n = 194).The black curve is the best-fit function for in situ AGB estimation, the black dashed line is the 45 • line, and the red solid line is the linear regression trend line.

Figure 5 .
Figure 5. Dynamic curves of rice dry aboveground biomass (g/m 2 ) based on logarithmic equations: (a) logarithmic curve of one sample plot in 2013; (b) comparison of measured aboveground biomass with estimated values of all the sample plots during 2012-2013 (n = 194).The black curve is the best-fit function for in situ AGB estimation, the black dashed line is the 45° line, and the red solid line is the linear regression trend line.

Figure 6 .
Figure 6.Regression analysis of HJ-1 CCD cumulative VIs versus dry aboveground biomass (g/m 2 ) using a quadratic polynomial fit function: (a) daily cumulative NDVI function and (b) 10-day composite data based on cumulative NDVI function.The blue solid lines are the regression trend lines.

Figure 6 .
Figure 6.Regression analysis of HJ-1 CCD cumulative VIs versus dry aboveground biomass (g/m 2 ) using a quadratic polynomial fit function: (a) daily cumulative NDVI function and (b) 10-day composite data based on cumulative NDVI function.The blue solid lines are the regression trend lines.

Figure 7 .
Figure 7. Results of dynamic mapping of LAI for single-cropped rice in Deqing County in 2012.Figure 7. Results of dynamic mapping of LAI for single-cropped rice in Deqing County in 2012.

Figure 7 .
Figure 7. Results of dynamic mapping of LAI for single-cropped rice in Deqing County in 2012.Figure 7. Results of dynamic mapping of LAI for single-cropped rice in Deqing County in 2012.

Figure 8 .
Figure 8. Results of dynamic mapping of LAI for single-cropped rice in Deqing County in 2013.

Figure 8 . 18 Figure 8 .
Figure 8. Results of dynamic mapping of LAI for single-cropped rice in Deqing County in 2013.

Figure 10 .
Figure 10.Dry aboveground biomass dynamic mapping of single-cropped rice in Deqing County in 2013.

Figure 10 .
Figure 10.Dry aboveground biomass dynamic mapping of single-cropped rice in Deqing County in 2013.

Figure 10 .
Figure 10.Dry aboveground biomass dynamic mapping of single-cropped rice in Deqing County in 2013.

Figure 11 .
Figure 11.Flow chart of rice growth dynamic mapping.

Figure 11 .
Figure 11.Flow chart of rice growth dynamic mapping.

Table 1 .
Dates of the field campaigns and corresponding HJ-1A/B charge-coupled device (CCD) images with sample numbers for each date.
Table 1．Dates of the field campaigns and corresponding HJ-1A/B charge-coupled device (CCD) images with sample numbers for each date.

Table 3 .
Correlation coefficients (r 2 ) between image features and rice growth parameters at different growth stages.

Table 4 .
Results of regression models at different single-cropped rice (SCR) growth stages.
E, P, and Q denote exponential, power, and quadratic polynomial fit of the traditional regression methods, respectively; B, S denote BPNN and SVM regression methods, respectively.