Estimation of rice phenology date using integrated HJ-1 CCD and Landsat-8 OLI vegetation indices time-series images

Estimation of rice phenology date using integrated HJ-1 CCD and Landsat-8 OLI vegetation indices time-series images 基于环境减灾卫星及Landsat-8卫星的植被指数时间序列的水稻物候期提取研究 1 first-author 1 2 2 1 1 1 1 1 1 corresp ‡Corresponding Author ‡Corresponding Author †E-mail:xdsongy@zju.edu.cn †E-mail:xdsongy@zju.edu.cn xdsongy@zju.edu.cn Institute of Remote Sensing and Information Application, College of Environmental and Resource Sciences, Zhejiang University, Hangzhou 310058, China Institute of Remote Sensing and Information Application, College of Environmental and Resource Sciences, Zhejiang University, Hangzhou 310058, China Institute of Remote Sensing and Earth Sciences, Hangzhou Normal University, Hangzhou 311121, China Institute of Remote Sensing and Earth Sciences, Hangzhou Normal University, Hangzhou 311121, China Accurate estimation of rice phenology is of critical importance for agricultural practices and studies. However, the accuracy of phenological parameters extracted by remote sensing data cannot be guaranteed because of the influence of climate, e.g. the monsoon season, and limited available remote sensing data. In this study, we integrate the data of HJ-1 CCD and Landsat-8 operational land imager (OLI) by using the ordinary least-squares (OLS), and construct higher temporal resolution vegetation indices (VIs) time-series data to extract the phenological parameters of single-cropped rice. Two widely used VIs, namely the normalized difference vegetation index (NDVI) and 2-band enhanced vegetation index (EVI2), were adopted to minimize the influence of environmental factors and the intrinsic difference between the two sensors. Savitzky-Golay (S-G) filters were applied to construct continuous VI profiles per pixel. The results showed that, compared with NDVI, EVI2 was more stable and comparable between the two sensors. Compared with the observed phenological data of the single-cropped rice, the integrated VI time-series had a relatively low root mean square error (RMSE), and EVI2 showed higher accuracy compared with NDVI. We also demonstrate the application of phenology extraction of the single-cropped rice in a spatial scale in the study area. While the work is of general value, it can also be extrapolated to other regions where qualified remote sensing data are the bottleneck but where complementary data are occasionally available. 目的 2 鉴于中国南方地区单季稻种植区在关键生育期较难获得清晰影像的情况，利用相互校准方法并结合HJ-1 CCD和Landsat-8 OLI传感器，生成时间分辨率更高的植被指数时间序列数据，并比较不同植被指数在提取水稻物候期中的差异。 创新点 2 本文通过传感器相互校准获得了具有更高时间分辨率的植被指数时间序列数据，同时研究了不同植被指数在提取水稻物候期中的有效性，从而提高了水稻物候期提取的精度。 方法 2 利用最小二乘法对HJ-1 CCD和Landsat-8 OLI传感器提取的植被指数（EVI2和NDVI）进行相互校准，验证了不同传感器可互补使用；利用一致性分析方法，对比不同植被指数在提取单季稻物候期中的有效性；通过极值法和最大斜率法提取研究区单季稻的移栽期、抽穗期和成熟期；将利用两传感器相结合形成的新植被指数时间序列数据得到的水稻物候期提取结果，与用单一传感器得到结果进行对比，分析水稻物候期提取精度。 结论 2 基于环境减灾卫星及Landsat-8卫星融合后得到的植被指数时间序列数据可以有效地提高南方单季稻物候期提取的精度，EVI2的提取效果优于NDVI，极值法和最大斜率法结合提取的单季稻物候期结果与野外调查及农气站统计结果较为吻合，可以很好地应用于实际业务中。 Phenological parameters Intercalibration Vegetation index HJ-1 CCD Landsat-8 OLI 物候提取 相互校准 植被指数 HJ-1 CCD Landsat-8 OLI Introduction 1 Dynamic variation of regional crop phenology is an important component of agriculture monitoring (Eerens et al., 2014 2014 ), and is of increasing relevance for environmental monitoring, e.g. changes in the phenological period and length of the growing season may be caused by climate variability (Brown and de Beurs, 2008 2008 ; Pan et al., 2015 2015 ). China has about one fifth of the world’s paddy rice land, and it is of critical importance to extract the phenology information of rice on a large spatial scale to serve related research in agricultural modeling, yield estimation, and climate change (FAOSTAT, 2012 2012 ; Li S.H. et al., 2014 2014 ). Remote sensing techniques offer a feasible tool for delineating spatiotemporal patterns of crop status on a per-pixel basis. Various kinds of optical remote sensing data and techniques have been applied in agricultural monitoring, e.g. the Advanced Very High Resolution Radiometer (AVHRR), the Moderate Resolution Imaging Spectroradiometer (MODIS), and the Systeme Probatoire d'Observation de la Terre (SPOT) VEGETATION (Wu et al., 2010 2010 ; Peng et al., 2011 2011 ; Cong et al., 2012 2012 ; Lanorte et al., 2014 2014 ). However, the spectral characteristics of the data may be seriously affected by a mixed-pixel problem due to the relative coarse resolutions (≥250 m) (Sakamoto et al., 2005 2005 ). The Landsat remote sensing data have relatively high spatial resolution but, hampered by its low-temporal resolution, are not suitable for phenological information extraction purposes in agriculture (Anderson et al., 2011 2011 ). When considering the influence of climate, the situation may be more serious. Most of the paddy fields located in the subtropical regions of China are influenced by the monsoon climate; as a consequence, the qualifying cloud-free remote sensing data during this monsoon period are relatively few (Pyongsop et al., 2010 2010 ; Cai et al., 2012 2012 ). To accurately extract the phenological information of rice in this region, a high spatial and temporal resolution remote sensing dataset is required, especially during the biologically sensitive growing period (Wu et al., 2010 2010 ). The small sun-synchronous satellites for environment and disaster monitoring and forecasting (HJ-1A and HJ-1B) from China were launched in 2008. The HJ-1A/B satellites have 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) as shown in Table 1 1 (She et al., 2015 2015 ). The sensor characteristics of HJ-1 CCD and Landsat-8 operational land imagers (OLIs), including the visible bands composition and spatial resolution, are very comparable. The Landsat-8 was launched on February 2013, and our intent is to combine HJ-1 CCD and Landsat-8 OLI data and verify its capability for extracting the phenological information of the single cropped rice in southern China. FOV: field of view; NIR: near infrared; Pan: panchromatic band Characteristic HJ-1 CCD Landsat-8 OLI  Altitude (km) 649 705  Repetition cycle (d) 31 16  Revisit interval (d) 4 (2 for constellation) 16  Scanning technique Push-broom Push-broom  Swath width (km) 360 (700 for two) 185  FOV (°) 0–35 15  Band ( n ) 4 9  Resolution (m) 30 15 (Pan), 30 3  Spectral interval (μm)  B1-Blue 0.43–0.52 0.45–0.51  B2-Green 0.52–0.60 0.53–0.59  B3-Red 0.63–0.69 0.64–0.67  B4-NIR 0.76–0.90 0.85–0.88 The time-series vegetation indices (VIs), such as the normalized difference vegetation index (NDVI) (Rouse et al., 1974 1974 ) and the 2-band enhanced vegetation index (EVI2) (Jiang et al., 2007 2007 ), are widely used in the studies of crop land classification, plant productivity, phenology, and crop growth monitoring (Panda et al., 2010 2010 ; Gao et al., 2013 2013 ; Shi et al., 2013 2013 ; Zhang et al., 2013 2013 ). It has been shown that VIs are relatively insensitive to the differences in angular viewing factors and atmospheric disturbances and thereby can be used as a benchmark for direct comparison between sensors (Steven et al., 2003 2003 ). Regression analysis has been common practice for the intercalibration of different remote sensing data. The ordinary least-squares (OLS) model is a primary tool for comparing two datasets and predicting one dataset from another (Ji et al., 2008 2008 ). Anderson et al. ( 2011 2011 ) compared the ResourceSat-1 NDVI values and Landsat-5 NDVI values by building an intercalibration equation. Li P. et al. ( 2014 2014 ) used OLS in cross-comparison of VIs derived from Landsat-8 OLI and Landsat-7 Enhanced Thematic Mapper Plus (ETM+) sensors, and showed that OLI and ETM+data are complementary. Agreement assessment is another important research topic in image intercalibration. The non-dimensional and symmetrical features of some agreement measures are more suitable for intercalibrating the VIs acquired from different sensors (Mielke and Berry, 2001 2001 ; Ji and Gallo, 2006 2006 ). In agreement analysis, the two datasets are treated equally and symmetrically, and the systematic and unsystematic differences between two datasets can be quantified. The geometric mean (GM) regression, orthogonal regression, and OLS bisector regression are widely used symmetrical regression models (Sprent and Dolby, 1980 1980 ; Isobe et al., 1990 1990 ; Rawlings et al., 1990 1990 ; Valsami and Macheras, 1995 1995 ). Ji et al. ( 2008 2008 ) used agreement analysis using symmetrical regression models and compared the AVHRR and MODIS 16-d composite NDVI products. GM regression-based agreement analysis was also used in a cross-sensor comparison study between Landsat-5 TM and IRS-P6 AWiFS (Chen et al., 2013 2013 ). The goal of this study was to generate a new VI time-series data with a higher temporal resolution by integrating HJ-1 CCD and Landsat-8 OLI images, in order to extract the phenological information of the single-cropped rice more accurately, and verify that VIs derived from HJ CCD and OLI images can be used as complementary data after proper intercalibration. We used the OLS method to integrate the two remote sensing datasets with the aid of a field campaign, and compared the efficiencies of NDVI and EVI2 using the agreement analysis technique. We also evaluated the integrated time-series VIs in rice phenological parameter extraction using field survey data. Data and methods 1 Study area 2 Deqing County located in North Zhejiang Province was selected as the study area. It has a mean annual temperature ranging from 13 to 16 °C and annual precipitation of 1379 mm (Fig. 1 1 ). This region is dominated by the tropical monsoon climate and more than 91% of the crop area is single-cropped rice, according to local statistical data. Two overlaid images are HJ-1B CCD1 (path 451 row 80) and Landsat-8 OLI (path 119 row 39), respectively. Dashed line refers to the HJ-1 CCD scene and the solid line refers to the Landsat-8 OLI scene (the swath width of each constellation) According to a field survey conducted during 2013–2014, the land cover types were classified as rice, trees, water bodies, artificial surfaces, and others. Four typical testing sites were selected with different land cover structures to intercalibrate the VIs derived from the two sensors. All testing sites were larger than 4 km×4 km. Furthermore, total 2281 independent randomly selected sample points were used for comparison based on the field survey and the Second National Soil Survey Vector Map, in which there were 484 for rice, 434 for trees, 709 for water bodies, 448 for artificial surfaces, and 206 for other areas. Sixteen site-based observations of rice phenology in 2013 provided by the Deqing agro-meteorological station were used as validation. Satellite data 2 The HJ-1 CCD sensors have similar characteristics to the Landsat-8 OLI, including a nominal descending orbit at 10:30 a.m. local crossing time, the same spatial resolution of 30 m, and overlapped bands in the visible spectrum (Fig. 1 1 ; Table 1 1 ). The prominent short revisit interval of HJ-1 CCD makes it an ideal candidate for crop status monitoring. We followed three basic principles in selecting remote sensing images of HJ-1 CCD and OLI in the study area to meet the requirement of calibration and time-series VI construction. They were that (1) the acquisition time difference of HJ-1 CCD and OLI should be within one day; (2) the time span of the remote sensing dataset should cover the whole growing season of the single-cropped rice; and (3) all the four testing sites should be cloud-free in the image pairs for calibration. Six pairs of archived HJ-1 CCD and Landsat-8 OLI data less than one day apart covering the study area in 2013 and 2014 were collected for intercalibration purposes (Table 2 2 ). To extract the rice phenological information, all the 32 relatively cloud-free HJ-1 CCD and OLI images which covered the rice growing period (from late May to early December) in 2013 were also collected to construct the time-series VIs (Table 3 3 ). GMT: Greenwich mean time. Sun azimuth of HJ ranges from 0° to 360°, while for Landsat-8 the range is from −180° to 180° Pair Sensor Acquisition date (year/month/day) Scene center scan time (GMT) Path/row Sun elevation (°) Sun azimuth (°) 1 HJ-1B CCD1 2013/04/14 02:06:46 448/80 58.01 302.94 Landsat-8 OLI 2013/04/14 02:33:23 119/39 61.02 131.11 2 HJ-1B CCD2 2013/07/18 01:46:37 448/80 59.57 279.88 Landsat-8 OLI 2013/07/19 02:33:33 119/39 67.17 108.62 3 HJ-1A CCD1 2013/08/05 02:07:42 451/80 60.86 289.65 Landsat-8 OLI 2013/08/04 02:33:35 119/39 65.27 116.26 4 HJ-1A CCD1 2013/11/09 01:56:45 450/80 34.85 323.72 Landsat-8 OLI 2013/11/08 02:33:19 119/39 40.30 157.66 5 HJ-1A CCD1 2013/12/10 01:56:09 451/80 29.44 326.70 Landsat-8 OLI 2013/12/10 02:33:10 119/39 33.61 157.88 6 HJ-1B CCD2 2014/01/27 01:39:26 451/80 30.06 319.56 Landsat-8 OLI 2014/01/27 02:32:41 119/39 35.73 150.42 No. Sensor Acquisition date (year/month/day) 1 HJ-1B CCD2 2013/05/14 2 HJ-1A CCD2 2013/05/20 3 HJ-1A CCD1 2013/05/28 4 HJ-1B CCD2 2013/06/06 5 Landsat-8 OLI 2013/06/17 6 Landsat-8 OLI 2013/07/03 7 HJ-1A CCD1 2013/07/05 8 HJ-1A CCD1 2013/07/09 9 HJ-1B CCD1 2013/07/11 10 HJ-1B CCD2 2013/07/18 11 Landsat-8 OLI 2013/07/19 12 HJ-1A CCD2 2013/07/24 13 HJ-1B CCD1 2013/07/30 14 Landsat-8 OLI 2013/08/04 15 HJ-1A CCD1 2013/08/05 16 HJ-1A CCD2 2013/08/12 17 HJ-1A CCD2 2013/08/16 18 HJ-1B CCD2 2013/08/18 19 Landsat-8 OLI 2013/08/20 20 HJ-1A CCD1 2013/09/01 21 HJ-1A CCD2 2013/09/04 22 HJ-1B CCD1 2013/09/26 23 HJ-1A CCD1 2013/10/11 24 HJ-1A CCD2 2013/10/20 25 Landsat-8 OLI 2013/10/23 26 HJ-1A CCD2 2013/10/28 27 Landsat-8 OLI 2013/11/08 28 HJ-1A CCD1 2013/11/09 29 HJ-1A CCD2 2013/11/16 30 HJ-1A CCD2 2013/11/20 31 HJ-1A CCD2 2013/11/28 32 HJ-1A CCD1 2013/12/02 The Level 2 HJ-1 CCD images were downloaded from the China Centre for Resources Satellite Data and Application (CRESDA). The Landsat-8 OLI standard level-one terrain-corrected (L1T) products were collected from the US Geological Survey (USGS). The OLI L1T products were geometrically corrected and orthorectified products by the data provider using standard systematic correction methods. Pro-processing and VI calculation 2 The images used were all processed by the radiometric calibration first. The formula and coefficients were collected from the raw data package. The atmospheric correction was performed using the Moderate Resolution Transmission (MODTRAN) 4 model integrated in FLAASH module in the ENVI package. The images were geometrically corrected using geographic information system (GIS) data to an accuracy less than 0.5 pixel (15 m). Two widely used VIs, i.e. NDVI and EVI2, were selected as the candidates to construct the time-series VI for extracting the phenological information of the single-cropped rice. The NDVI has been proved to be one of the best indicators for vegetation status monitoring (Xie et al., 2008 2008 ). While for EVI2, compared with EVI, it is more stable across sensors because only 2 bands are used (Jiang et al., 2007 2007 ; 2008 2008 ; Kim et al., 2010 2010 ; Wang et al., 2015 2015 ). The corresponding formulas are given as follows: $${\text{NDVI}} = \frac{{{\rho _{{\text{nir}}}} - {\rho _{{\text{red}}}}}}{{{\rho _{{\text{nir}}}} + {\rho _{{\text{red}}}}}}$$ , $${\rm{EVI2}} = 2.5 \times \frac{{{\rho _{{\rm{nir}}}} - {\rho _{{\rm{red}}}}}}{{{\rho _{{\rm{nir}}}} + 2.4 \times {\rho _{{\rm{red}}}} + 1}}$$ , where ρ red and ρ nir refer to the surface reflectance values of Bands 3 and 4 in HJ-1 CCD sensors, and Bands 4 and 5 in Landsat-8 OLI sensors, respectively (Table 1 1 ). Regression and agreement analysis 2 For intercalibration, we built the OLS regression functions and transformed Landsat-8 OLI-VIs to HJ-1 CCD-VIs based on the data of the four typical testing sites. The GM regression results were used in the following agreement analysis. The agreement metric, i.e. the agreement coefficient (AC), proposed by Ji and Gallo ( 2006 2006 ) was adopted to evaluate the agreement of the two datasets. The AC has the properties to distinguish between systematic and unsystematic differences. The AC is defined as: $${\text{AC}} = 1 - \frac{{\sum\limits_{i = 1}^n {{{({X_i} - {Y_i})}^2}} }}{{\sum\limits_{i = 1}^n {(|\overline X - \overline Y | + |{X_i} - \overline X |)(|\overline X - \overline Y | + |{Y_i} - \overline Y |)} }}$$ , where $$\bar X$$ and $$\bar Y$$ are the mean values of datasets X and Y , respectively. X i and Y i are the VI values of each pixel ( i =1, 2, 3, …, n ). The agreement between X and Y increases as AC approaches 1, and vice versa . The mean square difference (MSD), a metric used to measure the systematic and unsystematic differences between two dataset, is defined as: $${\text{MSD}} = \frac{1}{n}\sum\limits_{i = 1}^n {{{({X_i} - {Y_i})}^2}} {\text{. }}$$ The MSD can be further divided into the unsystematic mean product-difference (MPD u ) and systematic mean product-difference (MPD s ). The MPD u is as follows: $${\text{MP}}{{\text{D}}_{\text{u}}} = \frac{1}{n}\sum\limits_{i = 1}^n {(|{X_i} - {{\hat X}_i}|)(|{Y_i} - {{\hat Y}_i}|)}$$ , where $${\hat X_i}$$ and $${\hat Y_i}$$ are the predicted X and Y values, estimated by the GM and inverse GM regressions, respectively. The MPD s is then defined as follows: $${\text{MP}}{{\text{D}}_{\text{s}}} = {\text{MSD}} - {\text{MP}}{{\text{D}}_{\text{u}}}{\text{.}}$$ To compare the difference properties, i.e. the systematic and unsystematic differences, between the two datasets, we also calculated MPD s /MSD and MPD u /MSD, expressed as percentages. Rice phenological parameters extracted by time-series VIs 2 To evaluate the performance of time-series VIs produced by HJ-1 CCD and OLI sensors in rice phenological parameter extraction, we constructed two VI time-series for NDVI and EVI2, respectively. The 1st time-series VI was produced using HJ-1 CCD alone; and the 2nd one was generated by integrating HJ-1 CCD and OLI data. The Savitzky-Golay (S-G) filters were applied to smooth the VI time-series data to reduce the undesirable noise caused by unfavorable fluctuations (Jönsson and Eklundh, 2004 2004 ; Hird and McDermid, 2009 2009 ; Pan et al., 2015 2015 ). Readers interested in the S-G filters are referred to these studies (Jönsson and Eklundh, 2004 2004 ; Hird and McDermid, 2009 2009 ) for more details. The S-G filters were applied to the image-based VI time-series data composed by the HJ-1CCD and Landsat-8 OLI data from late May to early December in 2013 as listed in Table 3 3 . Before phenological parameters extraction, the planting area of single-cropped rice was estimated using a stepwise classification strategy proposed by Wang et al. ( 2015 2015 ). The total acreage of the single-cropped rice was about 94.0 km 2 in 2013, mainly distributed in the eastern part of the study area. Individual phenological stages of crops are identified by considering the VI increment or decrement between consecutive images over a certain period of time (Xin et al., 2002 2002 ; Boschetti et al., 2009 2009 ) (Fig. 2 2 ). Divided by the heading date, the growth phases of single-cropped rice can be classified into vegetative growth and reproductive growth. Because the maximum VI usually appears around the heading date, it is convenient to define the heading date as the date of the maximum VI on the VI profile. In general, the rice fields are flooded before transplanting and the VI of rice fields decreases during this period and then increases after rice planting. Therefore, it is reasonable to define the transplanting date of rice as the minimal point along the VI profile. Due to the etiolation and senescence of the rice leaves, the VI decreases after the heading, and the maturation date of rice is identified by the maximum slope method (Yu et al., 2003 2003 ; You et al., 2013 2013 ). The results were validated with the field survey data collected by the local agro-meteorological station in 2013. Smooth curve was generated using S-G filters Results 1 Intercalibration of VIs 2 We first compared the VI values of different land cover types of the sample points derived from HJ-1 CCD and Landsat-8 OLI images. It showed that EVI2 was more stable and comparable between the two sensors compared with NDVI (Fig. 3 3 ). This was further illustrated by the relatively small standard deviation of EVI2. The fluctuations of mean values of EVI2 were less than 0.04 while those of NDVI were around 0.05. dNDVI: differenced NDVI; dEVI2: differenced EVI2. dNDVI (a) and dEVI2 (b) are derived from Landsat-8 OLI and HJ-1 CCD images (image pair No. represents the six images pairs listed in Table 2 2 ) The linear regression relationships between Landsat-8 OLI-VIs and HJ-1 CCD-VIs for all samples (pixels) in the testing sites are shown in Fig. 4 4 . The two datasets are highly correlated for both NDVI and EVI2. The coefficients of determination of NDVI and EVI2 were 0.871 and 0.891, respectively. Black solid lines are the 45° lines and red solid lines are the linear regression trendlines (Note: for interpretation of the references to color in this figure legend, the reader is referred to the web version of this article) The agreement analysis showed that NDVI and EVI2 both had relatively high AC values (0.87 for NDVI and 0.89 for EVI2; Table 4 4 ). It’s notable that the MPD u values were much higher than the MPD s values for both VIs, indicating that the unsystematic difference was greater than the systematic difference. Additionally, the MPD u /MSD values were much higher than MPD s /MSD values for both VIs, an indication that the unsystematic difference was the primary difference between the two datasets. The MPD s and MPD u values of EVI2 were lower than those of NDVI, and it appeared that HJ-1 CCD and Landsat OLI may have more consistency in using EVI2. VI AC MSD (×10 −3 ) MPD s (×10 −3 ) MPD u (×10 −3 ) MPD s /MSD (%) MPD u /MSD (%) NDVI 0.87 7.11 0.77 6.34 10.82 89.18 EVI2 0.89 4.77 0.45 4.32 9.43 90.57 The systematic difference between HJ-1 CCDVIs and Landsat-8 OLI-VIs can be minimized using OLS. By taking Landsat-8 OLI-VIs as the dependent variables and HJ-1 CCD-VIs as the response, two linear regression functions for NDVI and EVI2 derived from the four testing sites within six image pairs could be used to rebuild the time-series data in the following phenological parameter extraction. Extraction of phenological parameters using integrated HJ-1 and Landsat-8 VI time-series images 2 We compared the smoothed NDVI and EVI2 time-series data at pixel level using S-G filters (Fig. 5 5 ). The phenology stages could be clearly identified from the VI profiles. Most of the noise was successfully eliminated from the VI time-series. The single-cropped rice usually is transplanted during late May to early July, and reaches the heading stage from late August to early September according to the weather and agricultural scheduling. By adding the Landsat-8 OLI data, the updated VI curves for both NDVI and EVI2 moved to the right. The transplanting, heading, and maturation dates of the single-cropped rice extracted by HJ-1 VI time series and the integrated HJ-1 and Landsat-8 VI time-series images in Deqing County are presented in Figs. 6 6 – 8 8 separately. For transplanting date, the dates extracted by NDVI time-series concentrated in late May to middle June; after integrating Landsat-8 OLI data, with more pixels moved into the 160–169 group. Similarly, the transplanting dates extracted by the EVI2 time-series data of HJ-1 concentrated in late May to early June, and the dates also delayed after integrating Landsat-8 OLI data. Since there were different rice varieties planted in Deqing, the heading dates were more dispersive from early August to early September, but the dates extracted by the NDVI time series showed an advancing trend compared with the EVI2 time-series. For maturation dates, both of NDVI and EVI2 showed less difference between the original dataset and the integrated dataset compared with other phenological stages. (a) HJ-CCD derived NDVI time-series; (b) Integrated NDVI time-series; (c) HJ-CCD derived EVI2 time-series; (d) Integrated EVI2 time-series (a) HJ-CCD derived NDVI time-series; (b) Integrated NDVI time-series; (c) HJ-CCD derived EVI2 time-series; (d) Integrated EVI2 time-series (a) HJ-CCD derived NDVI time-series; (b) Integrated NDVI time-series; (c) HJ-CCD derived EVI2 time-series; (d) Integrated EVI2 time-series The accuracy of the phenological dates extracted as above was evaluated using field-based phenology observations (Table 5 5 ; Fig. 9 9 ). The estimated dates for each phenological stage were generally within ±6 d compared with the observations. The differences between the two dates decreased after using the integrated VI time-series as source data for both NDVI and EVI2. Table 5 5 showed that the root mean square error (RMSE) of the estimated phenological dates using HJ-1 CCD images alone was larger than that of the integrated NDVI and EVI2 images. The RMSE of EVI2-derived phenology dates was smaller than that of NDVI in transplanting and heading dates. It showed that after integrating HJ-1 CCD and Landsat-8 OLI data, the VI time-series provided a more accurate estimation of rice transplanting, heading, and maturation dates than using single data sources. VI Data source Transplanting date (d) Heading date (d) Maturation date (d) 2 NDVI HJ-1 CCD 9.5 11.5 5.3 HJ-1 CCD+Landsat-8 OLI 5.7 10.4 4.5 2 EVI2 HJ-1 CCD 8.4 7.4 6.4 HJ-1 CCD+ Landsat-8 OLI 3.1 4.3 4.2 (a) Transplanting date extracted by NDVI time-series; (b) Transplanting date extracted by EVI2 time-series; (c) Heading date extracted by NDVI time-series; (d) Heading date extracted by EVI2 time-series; (e) Maturation date extracted by NDVI time-series; (f) Maturation date extracted by EVI2 time-series. DOY:day of year Discussion 1 Where remote sensing is used as a tool to facilitate the extraction of crop phenology information, "more data" is always the better policy. The crucial point in this practice is to combine the readily available remote sensing data, if possible, and generate a more dense time-series of data. In this study, we took full advantage of HJ-1 CCD, characterized by its high temporal resolution, and combined this with Landsat-8 OLI data to generate integrated time-series VIs and extract the key phenology information of the single-cropped rice. We selected NDVI and EVI2 in this study and compared their efficiencies in extracting the phenological information of rice. The comparison of VIs between HJ-1 CCD and Landsat-8 OLI showed slight differences due to the spectral band differences (Teillet and Ren, 2008 2008 ; Anderson et al., 2011 2011 ; Li P. et al., 2014 2014 ). As we had pointed out, the inconsistency between the two sensors could be attributed to systematic and unsystematic differences quantitatively described by agreement analysis. The systematic difference between the two sensors could be minimized by OLS, while the unsystematic difference could be effectively attenuated using VI, which could minimize the influence of the environmental factors. We built the empirical relationships between HJ-1 CCD-VIs and Landsat-8 OLI-VIs datasets at pixel level using OLS. By using the linear regression models, HJ-1 CCD-VIs and Landsat-8 OLI-VIs could be transformed interchangeably. Scatter plots of HJ-1 CCD-VIs versus Landsat-8 OLI-VIs for four testing sites within six image pairs demonstrated that the VIs obtained from HJ-1 CCD and Landsat-8 OLI were highly correlated, and EVI2 had a relatively high R2 compared with NDVI (Fig. 4 4 ). Not surprisingly, the agreement analysis also showed that EVI2 between the two sensors had significant lower systematic and unsystematic differences compared with NDVI. To extract the phenological parameters of rice, the S-G filters were applied to smooth the VI time series data, and we defined criteria to judge the dates of transplanting, heading, and maturation of the single-cropped rice by considering the VI increment or decrement between consecutive images over a certain period of time according to the field campaign work. Compared with the observed data, EVI2 showed a significantly lower RMSE than NDVI in the estimation of transplanting and heading dates, and this was testified by the relative higher agreement coefficient of EVI2 (Table 4 4 ). We noticed that NDVI showed the largest RMSE at heading date, and the heading dates estimated by NDVI advanced compared with the observations (Table 5 5 ; Fig. xref 9 9 ). Because NDVI is easily saturated in well-vegetated areas (Huete et al., 2002 2002 ), it can be seen that NDVI becomes insensitive at high values of leaf area index (LAI) (Gu et al., 2013 2013 ). Considering that the single-cropped rice usually reaches the peak of LAI between the ear differentiation and early heading stages, the advanced heading date estimated by using NDVI is not unexpected. Compared with NDVI, EVI2 is more robust in capturing the difference in well-vegetated areas and thus has more potential in phenology extraction. We combined extremum and maximum slope methods in rice phenology estimation. Due to the different methodologies applied in estimating rice phenology in field and remote sensing, an error between the two kinds of estimation is unavoidable. However, EVI2, especially the integrated EVI2, showed a significant consistency and agreement in phenology dates estimation. However, it should be noted that the empirical equation built between HJ-1 CCD-VI and Landsat-8 OLI-VI datasets was partially determined by the sample plots and dates of the image pairs. So, there are inevitable uncertainties in the estimation. If extrapolating this kind of integrating methodology and results to other regions or years and considering the variations of environmental factors and vegetation status, necessary re-calibration steps must be taken. Conclusions 1 Accurate extraction of the phenology information of rice in large spatial scale is crucial to agricultural management and related ecological studies. High quality time-series remote sensing data are of critical importance in identifying the key phenology dates as close as possible to the "true" dates. Most of the paddy fields in southern China are influenced by the monsoon climate, and qualified remote sensing data are scarce resources in these regions. In this study, we showed that by integrating HJ-1 CCD and Landsat-8 OLI data using OLS, the phenological parameters of the single-cropped rice can be estimated more accurately. Two widely used VIs, namely NDVI and EVI2, were adopted, and not surprisingly, the two indices obtained from HJ-1 and Landsat-8 showed high correlation and agreement (R2 \textgreater 0.87, AC \textless 0.86). However, compared with NDVI, EVI2 was more stable and comparable between the two sensors. Compared with the observed data, the integrated VI time-series had a relatively low RMSE, in which EVI2 was superior to NDVI. We also demonstrated the application of phenology extraction of the single-cropped rice in spatial scale in the study area. This work is of general value and can be extrapolated to other regions where qualified remote sensing data are the bottleneck and where complementary data are occasionally available. References JH Anderson , KT Weber , B Gokhale , 等 . Intercalibration and evaluation of ResourceSat-1 and Landsat-5 NDVI . Can J Remote Sens , 2011 . 37 ( 2 ): 213 - 219 . DOI: 10.5589/m11-032 http://doi.org/10.5589/m11-032 . Intercalibration and evaluation of ResourceSat-1 and Landsat-5 NDVI M Boschetti , D Stroppiana , PA Brivio , 等 . Multiear monitoring of rice crop phenology through time series analysis of MODIS images . Int J Remote Sens , 2009 . 30 ( 18 ): 4643 - 4662 . DOI: 10.1080/01431160802632249 http://doi.org/10.1080/01431160802632249 . Multiear monitoring of rice crop phenology through time series analysis of MODIS images ME Brown , KM de Beurs . Evaluation of multiensor semi-arid crop season parameters based on NDVI and rainfall . Remote Sens Environ , 2008 . 112 ( 5 ): 2261 - 2271 . DOI: 10.1016/j.rse.2007.10.008 http://doi.org/10.1016/j.rse.2007.10.008 . Evaluation of multiensor semi-arid crop season parameters based on NDVI and rainfall WW Cai , JL Song , JD Wang , 等 . High spatial-and temporal-resolution NDVI produced by the assimilation of MODIS and HJ-1 data . Can J Remote Sens , 2012 . 37 ( 6 ): 612 - 627 . DOI: 10.5589/m12-004 http://doi.org/10.5589/m12-004 . High spatial-and temporal-resolution NDVI produced by the assimilation of MODIS and HJ-1 data XX Chen , JE Vogelmann , G Chander , 等 . Cross-sensor comparisons between Landsat 5 TM and IRS-P6 AWiFS and disturbance detection using integrated Landsat and AWiFS time-series images . Int J Remote Sens , 2013 . 34 ( 7 ): 2432 - 2453 . DOI: 10.1080/01431161.2012.743690 http://doi.org/10.1080/01431161.2012.743690 . Cross-sensor comparisons between Landsat 5 TM and IRS-P6 AWiFS and disturbance detection using integrated Landsat and AWiFS time-series images N Cong , SL Piao , AP Chen , 等 . Spring vegetation green-up date in China inferred from SPOT NDVI data: a multiple model analysis . Agric Forest Meteorol , 2012 . 165 ( 0 ): 104 - 113 . DOI: 10.1016/j.agrformet.2012.06.009 http://doi.org/10.1016/j.agrformet.2012.06.009 . Spring vegetation green-up date in China inferred from SPOT NDVI data: a multiple model analysis H Eerens , D Haesen , F Rembold , 等 . Image time series processing for agriculture monitoring . Environ Modell Softw , 2014 . 53 154 - 162 . DOI: 10.1016/j.envsoft.2013.10.021 http://doi.org/10.1016/j.envsoft.2013.10.021 . Image time series processing for agriculture monitoring FAOSTAT . Statistical Database of the Food and Agricultural Organization of the United Nations . 2012 . Dec. 30, 2014 . Available from http://faostat3.fao.org/browse/Q/*/E http://faostat3.fao.org/browse/Q/*/E . Statistical Database of the Food and Agricultural Organization of the United Nations S Gao , Z Niu , N Huang , 等 . Estimating the Leaf Area Index, height and biomass of maize using HJ-1 and RADARSAT-2 . Int J Appl Earth Obs Geoinf , 2013 . 24 1 - 8 . DOI: 10.1016/j.jag.2013.02.002 http://doi.org/10.1016/j.jag.2013.02.002 . Estimating the Leaf Area Index, height and biomass of maize using HJ-1 and RADARSAT-2 Y Gu , BK Wylie , DM Howard , 等 . NDVI saturation adjustment: a new approach for improving cropland performance estimates in the Greater Platte River Basin, USA . Ecol Indic , 2013 . 30 ( 0 ): 1 - 6 . DOI: 10.1016/j.ecolind.2013.01.041 http://doi.org/10.1016/j.ecolind.2013.01.041 . NDVI saturation adjustment: a new approach for improving cropland performance estimates in the Greater Platte River Basin, USA JN Hird , GJ McDermid . Noise reduction of NDVI time series: an empirical comparison of selected techniques . Remote Sens Environ , 2009 . 113 ( 1 ): 248 - 258 . DOI: 10.1016/j.rse.2008.09.003 http://doi.org/10.1016/j.rse.2008.09.003 . Noise reduction of NDVI time series: an empirical comparison of selected techniques A Huete , K Didan , T Miura , 等 . Overview of the radiometric and biophysical performance of the MODIS vegetation indices . Remote Sens Environ , 2002 . 83 ( 1-2 ): 195 - 213 . DOI: 10.1016/S0034-4257(02)00096-2 http://doi.org/10.1016/S0034-4257(02)00096-2 . Overview of the radiometric and biophysical performance of the MODIS vegetation indices T Isobe , ED Feigelson , MG Akritas , 等 . Linear regression in astronomy . Astrophys J , 1990 . 364 ( 1 ): 104 - 113 . DOI: 10.1086/169390 http://doi.org/10.1086/169390 . Linear regression in astronomy L Ji , K Gallo . An agreement coefficient for image comparison . Photogramm Eng Remote Sens , 2006 . 72 ( 7 ): 823 - 833 . DOI: 10.14358/PERS.72.7.823 http://doi.org/10.14358/PERS.72.7.823 . An agreement coefficient for image comparison L Ji , K Gallo , JC Eidenshink , 等 . Agreement evaluation of AVHRR and MODIS 16-day composite NDVI data sets . Int J Remote Sens , 2008 . 29 ( 16 ): 4839 - 4861 . DOI: 10.1080/01431160801927194 http://doi.org/10.1080/01431160801927194 . Agreement evaluation of AVHRR and MODIS 16-day composite NDVI data sets ZY Jiang , AR Huete , Y Kim , 等 . 2-band enhanced vegetation index without a blue band and its application to AVHRR data . 2007 . Proc. SPIE 6679, Remote Sensing and Modeling of Ecosystems for Sustainability IV . San Diego, CA 667905 DOI: 10.1117/12.734933 http://doi.org/10.1117/12.734933 . 2-band enhanced vegetation index without a blue band and its application to AVHRR data ZY Jiang , AR Huete , K Didan , 等 . Development of a two-band enhanced vegetation index without a blue band . Remote Sens Environ , 2008 . 112 ( 10 ): 3833 - 3845 . DOI: 10.1016/j.rse.2008.06.006 http://doi.org/10.1016/j.rse.2008.06.006 . Development of a two-band enhanced vegetation index without a blue band P Jönsson , L Eklundh . TIMESAT–a program for analyzing time-series of satellite sensor data . Comput Geosci , 2004 . 30 ( 8 ): 833 - 845 . DOI: 10.1016/j.cageo.2004.05.006 http://doi.org/10.1016/j.cageo.2004.05.006 . TIMESAT–a program for analyzing time-series of satellite sensor data Y Kim , T Miura , Z Jiang , 等 . Spectral compatibility of vegetation indices across sensors: band decomposition analysis with Hyperion data . J Appl Remote Sens , 2010 . 4 ( 1 ): 043520 DOI: 10.1117/1.3400635 http://doi.org/10.1117/1.3400635 . Spectral compatibility of vegetation indices across sensors: band decomposition analysis with Hyperion data A Lanorte , R Lasaponara , M Lovallo , 等 . Fisher-Shannon information plane analysis of SPOT/VEGETATION Normalized Difference Vegetation Index (NDVI) time series to characterize vegetation recovery after fire disturbance . Int J Appl Earth Obs Geoinf , 2014 . 26 441 - 446 . DOI: 10.1016/j.jag.2013.05.008 http://doi.org/10.1016/j.jag.2013.05.008 . Fisher-Shannon information plane analysis of SPOT/VEGETATION Normalized Difference Vegetation Index (NDVI) time series to characterize vegetation recovery after fire disturbance P Li , LG Jiang , ZM Feng . Cross-comparison of vegetation indices derived from Landsat-7 Enhanced Thematic Mapper Plus (ETM+) and Landsat-8 Operational Land Imager (OLI) sensors . Remote Sens , 2014 . 6 ( 1 ): 310 - 329 . DOI: 10.3390/rs6010310 http://doi.org/10.3390/rs6010310 . Cross-comparison of vegetation indices derived from Landsat-7 Enhanced Thematic Mapper Plus (ETM+) and Landsat-8 Operational Land Imager (OLI) sensors SH Li , JT Xiao , P Ni , 等 . Monitoring paddy rice phenology using time series MODIS data over Jiangxi Province, China . Int J Agric Biol Eng , 2014 . 7 ( 6 ): 28 - 36 . DOI: 10.3965/j.ijabe.20140706.005 http://doi.org/10.3965/j.ijabe.20140706.005 . Monitoring paddy rice phenology using time series MODIS data over Jiangxi Province, China PWJ Mielke , KJ Berry . Permutation Methods: a Distance Function Approach . 2001 . New York 155 - 238 . DOI: 10.1007/978-1-4757-3449-2 http://doi.org/10.1007/978-1-4757-3449-2 . Permutation Methods: a Distance Function Approach Z Pan , J Huang , Q Zhou , 等 . Mapping crop phenology using NDVI time-series derived from HJ-1 A/B data . Int J Appl Earth Obs Geoinf , 2015 . 34 ( 0 ): 188 - 197 . DOI: 10.1016/j.jag.2014.08.011 http://doi.org/10.1016/j.jag.2014.08.011 . Mapping crop phenology using NDVI time-series derived from HJ-1 A/B data SS Panda , DP Ames , S Panigrahi . Application of vegetation indices for agricultural crop yield prediction using neural network techniques . Remote Sens , 2010 . 2 ( 3 ): 673 - 696 . DOI: 10.3390/rs2030673 http://doi.org/10.3390/rs2030673 . Application of vegetation indices for agricultural crop yield prediction using neural network techniques DL Peng , AR Huete , JF Huang , 等 . Detection and estimation of mixed paddy rice cropping patterns with MODIS data . Int J Appl Earth Obs Geoinf , 2011 . 13 ( 1 ): 13 - 23 . DOI: 10.1016/j.jag.2010.06.001 http://doi.org/10.1016/j.jag.2010.06.001 . Detection and estimation of mixed paddy rice cropping patterns with MODIS data RI Pyongsop , ZB Ma , QW Qi , 等 . Cloud and shadow removal from Landsat TM data . J Remote Sens , 2010 . 14 ( 3 ): 534 - 545 . . Cloud and shadow removal from Landsat TM data JO Rawlings , SG Pantula , DA Dickey . Applied regression analysis: a research tool . J Oper Res Soc , 1990 . 41 ( 8 ): 782 - 783 . DOI: 10.2307/2583482 http://doi.org/10.2307/2583482 . Applied regression analysis: a research tool JW Rouse , RH Haas , JA Schell , 等 . Monitoring the vernal advancement of retrogradation of natural vegetation . 1974 . Technical Report No. E74-10113, Texas AM University, College Station . TX, USA . Monitoring the vernal advancement of retrogradation of natural vegetation T Sakamoto , M Yokozawa , H Toritani , 等 . A crop phenology detection method using time-series MODIS data . Remote Sens Environ , 2005 . 96 ( 3-4 ): 366 - 374 . DOI: 10.1016/j.rse.2005.03.008 http://doi.org/10.1016/j.rse.2005.03.008 . A crop phenology detection method using time-series MODIS data B She , JF Huang , RF Guo , 等 . Assessing winter oilseed rape freeze injury based on Chinese HJ remote sensing data . J Zhejiang Univ-Sci B (Biomed & Biotechnol) , 2015 . 16 ( 2 ): 131 - 144 . DOI: 10.1631/jzus.B1400150 http://doi.org/10.1631/jzus.B1400150 . Assessing winter oilseed rape freeze injury based on Chinese HJ remote sensing data JJ Shi , JF Huang , F Zhang . Multi-year monitoring of paddy rice planting area in Northeast China using MODIS time series data . J Zhejiang Univ-Sci B (Biomed & Biotechnol) , 2013 . 14 ( 10 ): 934 - 946 . DOI: 10.1631/jzus.B1200352 http://doi.org/10.1631/jzus.B1200352 . Multi-year monitoring of paddy rice planting area in Northeast China using MODIS time series data P Sprent , GR Dolby . The geometric mean functional lationship . Biometrics , 1980 . 36 ( 3 ): 547 - 550 . DOI: 10.2307/2530224 http://doi.org/10.2307/2530224 . The geometric mean functional lationship MD Steven , TJ Malthus , F Baret , 等 . Intercalibration of vegetation indices from different sensor systems . Remote Sens Environ , 2003 . 88 ( 4 ): 412 - 422 . DOI: 10.1016/j.rse.2003.08.010 http://doi.org/10.1016/j.rse.2003.08.010 . Intercalibration of vegetation indices from different sensor systems PM Teillet , XM Ren . Spectral band difference effects on vegetation indices derived from multiple satellite sensor data . Can J Remote Sens , 2008 . 34 ( 3 ): 159 - 173 . DOI: 10.5589/m08-025 http://doi.org/10.5589/m08-025 . Spectral band difference effects on vegetation indices derived from multiple satellite sensor data G Valsami , P Macheras . The geometric mean functional relationship approach to linear regression in pharmaceutical studies: application to the estimation of binding parameters . Pharm Pharmacol Comm , 1995 . 1 ( 12 ): 551 - 554 . DOI: 10.1111/j.2042-7158.1995.tb00377.x http://doi.org/10.1111/j.2042-7158.1995.tb00377.x . The geometric mean functional relationship approach to linear regression in pharmaceutical studies: application to the estimation of binding parameters J Wang , J Huang , K Zhang , 等 . Rice fields mapping in fragmented area using multi-temporal HJ-1A/B CCD images . Remote Sens , 2015 . 7 ( 4 ): 3467 - 3488 . DOI: 10.3390/rs70403467 http://doi.org/10.3390/rs70403467 . Rice fields mapping in fragmented area using multi-temporal HJ-1A/B CCD images WB Wu , P Yang , HJ Tang , 等 . Characterizing spatial patterns of phenology in cropland of China based on remotely sensed data . Agric Sci China , 2010 . 9 ( 1 ): 101 - 112 . DOI: 10.1016/S1671-2927(09)60073-0 http://doi.org/10.1016/S1671-2927(09)60073-0 . Characterizing spatial patterns of phenology in cropland of China based on remotely sensed data YC Xie , ZY Sha , M Yu . Remote sensing imagery in vegetation mapping: a review . J Plant Ecol , 2008 . 1 ( 1 ): 9 - 23 . DOI: 10.1093/jpe/rtm005 http://doi.org/10.1093/jpe/rtm005 . Remote sensing imagery in vegetation mapping: a review J Xin , Z Yu , L van Leeuwen , 等 . Mapping crop key phenological stages in the north China plain using NOAA time series images . Int J Appl Earth ObsGeoinf , 2002 . 4 ( 2 ): 109 - 117 . DOI: 10.1016/S0303-2434(02)00007-7 http://doi.org/10.1016/S0303-2434(02)00007-7 . Mapping crop key phenological stages in the north China plain using NOAA time series images XZ You , JH Meng , M Zhang , 等 . Remote sensing based detection of crop phenology for agricultural zones in China using a new threshold method . Remote Sens , 2013 . 5 ( 7 ): 3190 - 3211 . DOI: 10.3390/rs5073190 http://doi.org/10.3390/rs5073190 . Remote sensing based detection of crop phenology for agricultural zones in China using a new threshold method FF Yu , KP Price , J Ellis , 等 . Response of seasonal vegetation development to climatic variations in eastern central Asia . Remote Sens Environ , 2003 . 87 ( 1 ): 42 - 54 . DOI: 10.1016/S0034-4257(03)00144-5 http://doi.org/10.1016/S0034-4257(03)00144-5 . Response of seasonal vegetation development to climatic variations in eastern central Asia LW Zhang , JF Huang , RF Guo , 等 . Spatioemporal reconstruction of air temperature maps and their application to estimate rice growing season heat accumulation using multi-temporal MODIS data . J Zhejiang Univ-Sci B (Biomed & Biotechnol) , 2013 . 14 ( 2 ): 144 - 161 . DOI: 10.1631/jzus.B1200169 http://doi.org/10.1631/jzus.B1200169 . Spatioemporal reconstruction of air temperature maps and their application to estimate rice growing season heat accumulation using multi-temporal MODIS data Compliance with ethics guidelines Jing WANG, Jing-feng HUANG, Xiu-zhen WANG, Meng-ting JIN, Zhen ZHOU, Qiao-ying GUO, Zhe-wen ZHAO, Wei-jiao HUANG, Yao ZHANG, and Xiao-dongSONG declare that they have no conflict of interest. This article does not contain any studies with human or animal subjects performed by any of the authors. 10.1631/jzus.B1500087 Project supported by the National High-Tech R & D Program (863) of China (No. 2012AA12A30703) and the Fundamental Research Funds for the Central Universities, China 2015-04-15 2015-09-17 2015-09-14 2022-01-05 Copyright © Zhejiang University and Springer-Verlag Berlin Heidelberg 2015 2015 Journal of Zhejiang University-SCIENCE B (Biomedicine and Biotechnology) 10 16