A Study on Phenology Detection of Corn in Northeastern China with Fused Remote Sensing Data A Study on Phenology Detection of Corn in Northeastern China with Fused Remote Sensing Data

Accurate phenology information detection is the basis for other remote-sensing based agriculture applications. So far, there have been a lot of phenology estimation models based on remote-sensing data, but little attention was paid to microscopic mechanism of crops and the environmental factors. The main purpose of this chapter is to apply a new phenology detection model, which combined physical mechanism-based crop models with remote-sensing data to detect the critical phenological stages of corn in Northeast China (Jilin and Liaoning Provinces). Compared to the phenology observations from the agriculture meteorological stations, the corn phenology estimation accuracy in Northeast China using only MODIS data is much lower than that in the US field sites. The main reason might be the small size of single piece of cropland in northeastern China, which led to the mixed MODIS pixels. Accordingly, Landsat and MODIS data fusion methods were applied to get time-series images with Landsat-like spatial resolution and MODIS-like temporal resolution, and quantitative and qualitative validation was conducted to evaluate and verify the accuracy of the data fusion. The results show that data fusion of Landsat and MODIS improved the spatial resolution and decreased the inﬂuence of mixed pixels.


Introduction
Accurate measurements of regional-to-global-scale vegetation dynamics and phenology information improve our understanding of inter-annual vegetation change in terrestrial ecosystems, as well as climatic and other environmental variations from year to year [1][2][3][4][5][6][7]. The phenological stages of crops provide essential information for agricultural activities such as irrigation scheduling and fertilizer management [26]. In addition, accurate detection of key crop growth stages is a key input for crop yield estimation based on remote-sensed vegetation index (VI) data and crop model applications integrated with remote sensing data [8][9][10][11].
The traditional approaches to estimating crop phenology have been through ground observations and the use of crop models (e.g., Simple and Universal CROp growth Simulator (SUCROS) [12], Hybrid-Corn [14], WOrld FOod STudies (WOFOST) [13][14][15]). The crop models can estimate crop growth dates with a high level of accuracy [root mean square error (RMSE): 0-4 days], but they require a number of detailed information inputs such as crop (e.g., cultivar used and plant population), weather (e.g., temperature, rainfall, solar radiation and wind speed) conditions [14,16] and soil (e.g., initial soil moisture). On the one hand, the use of these crop models is usually limited by the availability of the required data inputs; on the other hand, the models need to be calibrated for particular species and site-specific conditions based on ground data [11,16], and the ground observations collected by observers are not cost-efficient. Accordingly, the traditional methods are sitespecific and typically cannot monitor crop phenology beyond the field scale over larger areas [11].
Satellite remote sensing observations from global imaging sensors offer considerable potential to provide the information of regional spatio-temporal patterns of the key crop growth stages of cash crops in a consistent, time-and cost-efficient manner [11]. The commonly used remote-sensing based phenology detection methods can be divided into four groups [11]: (1) threshold methods that estimate phenological stages by using either a fixed or dynamic threshold [17][18][19][20]; (2) moving window methods that determine the phenology dates by vegetation index (VI) changes of a time-series VI curve in a defined moving temporal window (e.g., 20 days) [21][22][23]; (3) function fitting methods that apply mathematical functions (e.g., logistic [24], Fourier transformation, wavelet [25]) to fit the time-series VI curves to a given function and extract phenological stages through the detection of defined feature points (e.g., second derivative equals 0) on the function curves; and (4) the shape model method [two-step filtering (TSF) approach developed by Sakamoto et al. [26]] that applies a novel shape-model fitting concept to times-series VI curves for date identification and by Zeng et al. [11] who proposed a hybrid method with environmental factors taken into consideration by integrating crop models.
The first three remote-sensing-based phenology detection methods summarized above are generally based on mathematical methods that directly detect the feature points as the transition dates of vegetation [24,27,28]. Usually, these dates represent general vegetation growth stages (e.g., greenup onset, peak greenness and dormancy onset [24]) but have little association with the specific agronomic stages of a specific kind of crops (e.g., corn, soybeans) [11]. Some key crop growth stages (e.g., R1 (beginning bloom), R3 (beginning pod) stages of soybeans) are challenging or impossible to be detected by finding feature points of the VI time-series curve. In addition, these methods are often sensitive to observation errors and noise caused by cloud coverage, atmospheric constituents (e.g., water vapor), bi-directional reflectance distribution function (BRDF) effect and the mixed-pixel effect in time-series VI data products [11,26].
The TSF approach is potential for the specific agronomic stages detection of specific crops with a high estimation accuracy [root mean square error (RMSE) ranging from 2.9 to 7.0 days and from 3.2 to 6.9 days for corn and soybeans, respectively] [26]. However, this method was based on an assumption that the shape model is linearly scalable to fit the time-series VI curve through geometrical scaling, regardless of all the factors that would influence the crop's growth pattern expressed in the VI data [11,26]. For example, air temperature is generally one of the decisive factors that affect the growth rate of all crops. Photoperiod is even a more important factor for some photosensitive crops, as longer daylength may decrease the development rate by delaying reproductive development (e.g., soybeans). As a result, the time-series VI profile of a crop's growth pattern under different environmental conditions can vary from year to year [11].
Zeng et al. [11] proposed a hybrid phenology detection approach that incorporates the "shapemodel fitting" concept of the TSF method [26] and simulation concept of traditional crop model that incorporates other environmental factors that influence crop development [11]. The approach is designed to detect the critical vegetative stages (V1 and V6) and reproductive stages (R1-R6) of corn from the MODIS 250-m Wide Dynamic Range Vegetation Index (WDRVI) time-series data. This method was tested over a 10-year period (2003-2012) for three experimental field sites to calibrate and quantitatively assess its performance with groundbased crop phenology observations for each site. It was also tested regionally over eastern Nebraska and the state of Iowa to evaluate its ability to characterize spatio-temporal variation of the targeted corn and soybean phenology stage dates across a larger major crop-producing region.
In this chapter, this method was tested in northeastern China where the environment is different with that in the US. As the piece of farmland is smaller than that in the US, the results were affected by mixed pixels from the MODIS observations. Then data fusing method was used to generate time-series VI data with high spatial and temporal resolution data, based on which phenological information of corn was estimated.

Study area
Northeast China, an important agricultural production base, plays a significant role in the national food security. In addition, it is a major corn producing area and named as one of the world's three golden Corn Belts. Jilin and Liaoning province were selected as the study area to further study the phenology estimation in nonirrigated region with small plots of cropland.
The ground-based data are from 30 agricultural meteorological stations (Figure 1). It mainly includes crop species, growth stage, soil moisture and air temperature. The ground-based crop growth stage observations were conducted by agronomists once every 1-10 days during the growing seasons and recorded in every 10 days. The recorded stages include VE, VT, R1, R4, R6 stages if they were match to general stage system (vegetation (V)-stage or reproductive (R)-stage) according to their growth stage names.

WDRVI data
The original source data for the WDRVI time-series data set used in this study include the MODIS 8-day composite, 250 m surface reflectance data (MOD09Q1, Collection 5) and 30 m Landsat Surface Reflectance Climate Data Records (CDRs). The WDRVI developed by Gitelson [29] offers valuable alternative VI to monitor crops because its sensitivity is at least three times greater than that of NDVI at moderate-to-high LAI values and maintains a linear relationship across the range of LAI values for both corn and soybeans [11,26,29,30].

Average air temperature data calculated from MODIS land surface temperature products
Average air temperature is a critical input for crop models in this chapter. It is estimated by both daily daytime and nighttime 1-km land surface temperature (LST) data (MOD11A1, Collection 5) by a linear regression model [11,31], and the calibration and validation were based on the observations from the meteorological stations in the study areas of China. The daily averaged air temperature (Tavg) was calculated by averaging the daily maximum and minimum air temperatures [11]. The 8-day composite air temperature was then calculated by averaging daily Tavg of cloud-free days during the 8 day composite period [11].

Land use data layer
GlobeLand30 datasets were selected as the regional land use classification data of China, which was developed under the "Global Land Cover Mapping at Finer Resolution" project led by National Geomatics Center of China (NGCC) China [32]. An area-ratio threshold of 75% was adopted to select all MODIS 250 m pixels across the region that was covered predominately by corn [11,26].

Using spatial and temporal adaptive reflectance fusion model algorithm to blend Landsat and MODIS surface reflectance
Before data fusion, preprocessing of Landsat and MODIS data was conducted, including format conversion, atmospheric correction, geometrical correction, resampling, reprojection, registration, mosaicking and cropping. Using more than one pair of reference images, particularly those observed in different stages of the growing season, can reduce the impact of the cloud and improve the estimation accuracy [33]. Using the reference images in different stages of the growing season may avoid some mistakes caused by harvest and fallow. In addition, using the reference images that were observed closer to the time of the Landsat images to be estimated can improve the estimation accuracy [34]. In this chapter, three cloud-free Landsat images observed in the beginning, middle and end of growing season respectively were selected as the reference images as possible to be blended with 250 m 8-day MODIS reflectance image using spatial and temporal adaptive reflectance fusion model (STARFM) method. The two reference images observed in the beginning and middle of the growing season were used to estimate time-series high spatial resolution images of the beginning and middle of the growing season. Correspondingly, the two reference images observed in the middle and end of the growing season were used to estimate timeseries high spatial resolution images of the middle and end of the growing season. The quality flag data and cloud mask of both MODIS and Landsat were used to identify the pixels with cloud or low quality.

Crop models for corn
Plant growth processes are mainly influenced by interactions among genotype, environment conditions and crop management [35,36]. In this chapter, we assumed that the influence of genotype and crop management (i.e., proper management of pests and diseases and fertilizer applications was implemented) was minimal as compared to environmental factors [11]. For this chapter, the developmental stages of corn were assumed to be most closely related to air temperature [37][38][39].
The models that describe the relationship between crop development and environment factors include linear and nonlinear models. Nonlinear approaches have been shown to provide better predictions of plant development stages than linear models for a number of different crops including corn [40], soybeans [15], wheat [41], potato [42], rice [43], as nonlinear models often describe the biological processes underlying crop growth in more detail [44,45].
The Wang-Engel (WE) model [45], nonlinear models, simulates crop development with response functions that range from 0 to 1. The temperature response function in the WE model is described by a beta function with three parameters (i.e., minimum, optimum and maximum air temperatures). When the temperature is below the minimum or above the maximum temperature, crop development stops and the temperature response function equals zero. When the temperature is at optimum, which is a value between minimum and maximum temperature (e.g., 28°C for corn), development takes place at the maximum rate and temperature response function equals one [45]. The WE model was originally developed for winter wheat [45] but has also been used to simulate the development of other annual crops including corn [44] and soybeans [15] with positive results. Accordingly, the WE model was used in this study to describe the development of corn responding to air temperature Eqs. (1)-(4). (1) where r max is the maximum development rate (per day). T is the average near surface air temperature estimated from MODIS data. f(T) is a temperature response function, which varies from 0 to 1.T up , T base and T opt are the abovementioned three parameters (maximum, minimum, and optimum air temperature, respectively).

Combining crop models with time-series WDRVI data
Photothermal time which combined both temperature and photoperiod information was used to describe leaf appearance rate and phenological response of various plant species [46][47][48].
To combine the crop model and time-series MODIS WDRVI data, photothermal time (accumulated photothermal time, APTT) was used instead of calendar time [Day of Year (DOY)] on the time axis in this chapter. It was defined as the accumulated development rate (∑r) calculated by temperature and photoperiod response function (Eq. 5). Usually, the planting date varied from year to year to ensure that APTT value was calculated from the same onset of growth, a starting point of APTT was set as the beginning of greenness of the crops (S 0 ) derived by the Zhang's method [24].
where onset is when corn reaches the beginning of greenness detected by the Zhang's method [24]. f(T) and f(P) are temperature and photoperiod response functions, respectively. While as the growth of corn is insensitive to photoperiod to some extent, f(P) equals constant 1. In order to combine APTT time and time-series WDRVI data, the APTT was calculated at the same eight-day interval with the eight-day composite periods of the MODIS WDRVI data [11].

Building a shape model
Before building the shape model, the Savitzky-Golay Filter method [49] was used to de-noise and rebuild the time-series VI data. In order to generate the shape model from an idealized temporal VI for both crops, the shape model in this study built under APTT with higher peak WDRVI value (>80). In order to generate the shape model from an idealized temporal VI for corn, the smoothed WDRVI time series data were stretched in the Y direction to make sure that the peak WDRVI value equals the median peak value and shifted in the X direction by detecting the beginning of the growing season with the Zhang's method to make sure that the time-series curves have the same onset. After the preliminary Y-scale and X-shift, the discrete points of smoothed WDRVI data were used to build the shape models for corn. They were fitted by the sum of three sine functions [11]. The predefined APTT X 0 for each phenological stage was calculated by averaging the APTT of each phenological stage' transition dates of all the observations used to build the shape model [11].

Fitting the smoothed WDRVI data on the shape model
The smoothed WDRVI time series data for each year were stretched by a Y-scale to make sure that the peak value was equal to the median peak value. Then, the shape model is fitted on the smoothed WDRVI time-series data. The optimum scaling parameters (X-scale, Y-scale) that approximate the fit of the smoothed WDRVI data to the shape model were calculated based on the smallest root mean square error (RMSE) between the shape model and the scaled smoothed WDRVI data [11], see Eqs. (6) and (7).
The function h(x) refers to the smoothed WDRVI data, f(x) is the shape model and g(x) is transformed from the smoothed WDRVI data for a given site or year, where x is the APTT value.

Shape model
The smoothed WDRVI curve with peak larger than 80 was selected to establish the corn shape model in the study area to reduce the impact of mixed pixels, noises, environmental stress on crop growth and phenology (Figure 2).

Phenology estimation results
Due to the lack of regional-scale crop phenological observations in Northeast China, phenology estimation accuracy was only validated in a field scale by comparing with the observations from agricultural meteorological stations and the estimated phenological stages derived from the time-series WDRVI data of the corn field in the stations. The RMSE of the estimated five corn growth stages varies from 3.78 to 8.41 days with little system bias (Figure 3 and Table 1). The estimation accuracy was higher in the field sites of University of Nebraska (with RMSE varied from 1.99 to 4.30 days) [11]. The RMSE of the estimated five phenological stages is within 15 days, which indicates that the phenology estimation model is also effective in northeastern China.
Similar to the results in the field sites of the University of Nebraska Lincoln, lower estimation accuracy was observed at the beginning and end of the growing season (such as V1, R6) than that at the middle stage of the growing season when vegetation cover reaches the peak, such as R1 stage ( Table 1). As in the beginning and end of the growing season, the vegetation index is more sensitive to noise and susceptible to weeds and other noncrop plants [50,51].

The comparison between the model from northeastern China and the field sites from Nebraska
The shape models of corn which were built on the data from the field sites of University of Nebraska-Lincoln and the agricultural meteorological stations in northeastern China were compared and shown in Figure 4. The needed APTT time of corn in two study areas was about the same, while there were some differences between the derived time-series VI curves from the study areas in the USA and China [11]. The time-series VI curves from the field sites of University of Nebraska-Lincoln and the Corn Belt in the USA were generally in agreement with each other: the curves reach the peak and then drop quickly after a slow decline, while the curves derived from corn field in northeastern China drop rapidly after reaching the peak [11].  The diameters of the farmland in the Corn Belt of USA are usually greater than 250 m, the size of MODIS pixels in red and near-infrared bands, and the VI observations from MODIS were less influenced by mixed pixels. Therefore, the time-series VI curves from this area may be closer to reality, which decrease rapidly after a slow decline following the peak. While in northeastern China, the size of the cropland is usually smaller than 250 m, the mixed pixels may cause error to time-series VI curves. Although the shape of the time-series VI curves in some stations of northeastern China is to some extent consistent with those in the Corn Belt of USA, the shape model from northeastern China was built by fitting the data with mixed pixels and pure pixels, which results in the loss of some characteristic details of the time-series VI curves.
The predefined APTT time required for some phenological stages in the field of the USA and China is about the same, such as the R1 and R6 stages, while there are obvious differences at some stage. The needed APTT time for the corn from the Corn Belt of USA reached V1 stage even less than the needed APTT time for VE stage of corn from the agricultural meteorological stations in northeastern China. One reason may be the description of phenology for corn in the two regions is different. In addition, the differences in the corn varieties and genotype may also contribute to difference of the need APTT.

Shape model
As in northeastern China, the size of the cropland is usually smaller than 250 m, the mixed pixels may cause error in time-series VI curves. Data fusing of 8-day composition MODIS images and Landsat images was conducted with the StarFM method to get the 8-day composition data with the spatial resolution like Landsat and with the temporal resolution like MODIS. The smoothed WDRVI curve of corn field in each meteorological station was calculated from time-series fused data and MODIS data, respectively. The smoothed WDRVI curves with peak larger than 80 in 2013 were selected to establish the corn shape model in the study area to reduce the impact of mixed pixels, noises, environmental stress on crop growth and phenology ( Figure 5).

Phenology estimation results
The error of five corn growth stages estimated from fused data is basically within 10 days (Figure 6). The RMSE of the results varies from 2.63 to 6.23. The estimation accuracy is higher than that derived from MODIS data only (3.78-8.41 days) ( Table 2). Similar to the results in the field sites of the University of Nebraska Lincoln [11], lower estimation accuracy was observed at the beginning and end of the growing season (such as V1, R6) than that at the middle stage of the growing season when vegetation cover reaches the peak, such as R1 stage ( Table 2).

The comparison between the model from the MODIS only data and the data fusion of Landsat and MODIS
The estimation accuracy of five corn growth stages from MODIS data only is lower than that derived from MODIS and Landsat fused data. As the mixed pixels in MODIS data may include artificial surface, water bodies, fallow land, etc. around the cropland, which may result in a decrease in the WDRVI data of the entire MODIS pixel where the station locates and influence the phenology estimation accuracy. Similarly, the forests area around the cropland might increase the WDRVI value in the early and late stages of crop growth. By sampling 9 × 9 window of the 30 m pixels (270 × 270 m) centered on the farmland of the station, Figure 7 shows the time-series WDRVI curves derived from MODIS data only and the fused data in 2013. The trend of the fused time-series WDRVI data is the same with that of MODIS derived data, as cloud and other noise affect the MODIS data as well as the fused result with MODIS data and Landsat data. However, the fusion of MODIS and Landsat data can reduce the influence of mixed pixels to a certain extent. For example, station 1 is mixed with cultivated land and artificial surface, in which the proportion of artificial surface in the two sites 9 × 9 windows was 50.6%. The fused data of this site are higher than the observed values of MODIS images.   Station 2 and station 3 are mixed with cultivated land and fallow land. The fused data of these two sites are higher than that obtained from MODIS images, especially in the middle of growing stage when WDRVI value is high in cropland, but low in fallow land.

Conclusion
This chapter applied a hybrid remote sensing-based crop phenology estimation method for corn that incorporated the simulation concept of crop growth models with the shapemodel fitting concept of the TSF method developed by Sakamoto et al. [26], in Northeast China (Jilin and Liaoning Provinces), which is an important grain production base in China. Compared to the field phenology observations from the agriculture meteorological stations, the corn phenology estimation accuracy in Northeast China using only MODIS data is much lower than that in the US field sites (RMSE of corn phenology estimate in Northeast China ranges from 3.78 to 8.41 days). The main reason might be the small size of single piece of cropland in northeastern China, which led to the mixed MODIS pixels. Accordingly, Landsat and MODIS data fusion methods were applied to get time-series images with Landsat-like spatial resolution and MODIS-like temporal resolution, and quantitative and qualitative validation was conducted to evaluate and verify the accuracy of the data fusion.
The results show that data fusion of Landsat and MODIS ensured the temporal resolution of time-series images, to some extent, improved the spatial resolution and decreased the influence of mixed pixels.