IntroductionTop

Wheat (Triticum aestivum L.) is one of the most important cereals grown in the world. More humans consume wheat as their main food more than other cereal grains (Pimentel & Pimentel, 2007). So, wheat plays a very important role in the world’s supply chain and food security. Therefore, maximum yield achievement at minimum expense is the main goal of wheat production. Awareness of crop growing dates, such as planting and harvesting dates, helps farm managers to reach the objectives and farm management. Harvest date is an important factor in field management. Planning the crop appropriate harvest date helps to schedule for harvest operation with less cost and maximize the profit (Abawi, 1993; Suwannachatkul et al., 2014). In these studies, wheat appropriate harvest date (WAHD) is the date in which the wheat has reached maturity, and its moisture content is convenient for harvest by a combine harvester. In Iran, the appropriate moisture content is in the range of 12-14 percent (Mansouri-Rad, 2000). If wheat is harvested earlier than WAHD, there will be yield loss and its moisture content will be inappropriate for harvesting and storage (Burnett & Bakke, 1930; Philips & O'Callaghan, 1974; Sabir et al., 2005). In contrast, harvesting after this date causes dry matter yield diminution and grain quality degradation. Moreover, yield loss due to weather conditions (such as temperature, rain, and wind), birds and insects damages to the matured product are disadvantages of delay in the harvest (BoIIand, 1984; Abawi, 1993; Farrer et al., 2006; Sun et al., 2007). In some areas of the world, the majority of farmers are smallholders, and the size of most agricultural farms is small. Consequently, it is not affordable for them to own their combine harvesters, leading to a shortage of harvesters in the harvest time (Bougari et al., 2013). Lack of necessary machinery during harvest campaign could further delay the harvest date. Therefore, the wheat producers require prior planning for providing the combine harvester rental and labors which requires former notice of WAHD for each field. Field surveying operations are costly, time-consuming and infeasible in large areas to determine harvest date (Moran & Pearce, 1997; Pinter et al., 2003). Some researchers (Porter & Gawith, 1999; McMaster& Wilhelm, 2003; Streck et al., 2003; Evers et al., 2010) have improved wheat phenological stages prediction models using photoperiod, water and vernalization to predict wheat maturity date. Other models such as STICS (Simulateur mulltidiscplinaire pour les Cultures Standard) and WOFOST (Worlds Food Study) have been developed by taking more factors into account such as temperature, nutrient and water stress for crop growth stage prediction (Boogaard et al., 1998; Brisson et al., 1998). Most of these models are based on meteorological parameters, making it difficult to discriminate within and between field differences in a dense configuration of the spatial grid (Meng et al., 2015). Moreover, the main limitation of crop models is the challenges in preparing reliable input data. The uncertainty about the spatial repartition of soil properties and micro meteorological variables at farm scale limits these model output’s assurance (Therond et al., 2011).

Remote sensing provides appropriate and timely images of the agricultural farms. The higher revisit frequency capability is the merit of remote sensing (Atzberger, 2013) for collecting the farm information. In recent decades, satellite remote sensing has been applied for agricultural field management operations such as yield and biomass estimation (Panda et al., 2010; Ren et al., 2008; Xie et al., 2009), phenological date prediction (Sakamoto et al., 2010; De Bernardis et al., 2016), and mapping of land use (Galford et al., 2008; Atzberger & Rembold, 2013). Vegetation indices (VIs), which are computational combinations of different spectral bands of the electromagnetic spectrum, simplify the analysis and processing of big data obtained by satellites (Govaerts et al., 1999; Viña et al., 2011). The strong contrast of absorption and scattering of the red and near-infrared bands can be combined into different quantitative indices to explain the vegetation conditions. VIs are the semi-analytical measurements of plant vegetation activity which ex­plain vegetation conditions during the growth stages. The benefit of VIs utilization is spectral reflectance data enhancement, by considering the variability of vegetation and minimizing of the atmospheric effect, sun-target-sensor geometry and soil (Moulin, 1999; Viña et al., 2011).

The relationship between remote sensing indices and plant bio-physical variables is nonlinear (Haboudane et al., 2004). Therefore, the usage of methods creating a non-linear relation between independent variables (VIs) and dependent variables could help predict better the plant bio-physical behavior changes. The artificial neural network (ANN) is the technique, widely used in the field of geo/bio-physical variable detections (Beale et al., 2008). Application of ANN techniques using the VIs and also visible blue, green, red, near-Infrared (NIR) and short wave (SWIR) regions of the electromagnetic spectrum have led to successful results for crop monitoring, crop cover, crop growth, crop nitrogen, crop yield and biomass estimation (Chen & McNairn, 2006; Karimi et al., 2006; Li et al., 2007; Xie et al., 2009; Prasad et al., 2012). The aim of this study was to generate a harvest schedule regional (HSR) map for predicting WAHD by using two years (2015 & 2016) VIs data derived by operational land imager (OLI) sensor (Landsat-8) and ANN techniques.

Material and methodsTop

Study site

The study was carried out in Qorveh plain in sou­theast of Kordestan province, West of Iran (35°15´N, 47°80´E), altitude 1900 m, with cold snowy winters and temperate summers (Fig. 1). The mean annual temperature of the site is 10.6 °C, and the average annual precipitation is 439 mm. The rainfed winter wheat is the major crop in this region, and it is sown in mid-October and harvested in mid-July of next year. The wheat-growing season is about 250 days.

Figure 1. Study site and sampling farm locations. Background: Landsat-8 OLI at 17 May 2016 (Row: 36, Path: 166).

Data acquisition

Field data

Field data for rainfed wheat yield measurement was collected during the wheat maturing stage of 2015 and 2016 years (from June 19 to July 11). Fifty sampling plots were selected from fifty wheat farms, in each year. They were flat and homogenous with areas bigger than 175×175 m2. The size of plots chosen was 60×60 m2 and plots were located at the center of wheat farms. A GPS receiver (Garmin GPSMAP 62s, Taiwan, with a spatial accuracy 3 m) was used for recording the geo-coordinates of each plot. The plant density in each plot was obtained in early June using a 1 m × 1 m quadrate. The plant density measurements were done by ten times random throwing of the quadrate through the zigzag sweep path and counting the number of plants in each throw and averaging of plant numbers of the quadrates. According to the field observations, plant densities often were constant on each farm.

The yield sampling was carried out by five random throws of the quadrate and choosing three plants in each throw. The plants were randomly chosen in the quadrate. At each sampling, fifteen plants were sampled for each plot, in total. By measuring the average grain weight of sampled plants in each plot and knowing plant density, the yield amount was acquired for each plot. The grains were used to measure the yield of fully matured plants. The sampling operation was conducted at a frequency of 2 days, from June 19 to July 11. Generally, eleven yield sampling observations were carried out for each plot in each year. The interpolation between yield observations was done by fitting a spline curve. The WAHD was considered the day of wheat highest yield.

Satellite image acquisition and pre-processing

Landsat-8 is the latest in a series of Landsat satellites, and it was launched on February 11, 2013 (http://science.nasa.gov/missions/ldcm). The operational land imager (OLI) and the thermal infrared sensor (TIRS) are two sensors which are carried by the Landsat-8 (Roy et al., 2014). We used OLI sensor data in this study. The appropriate spatial resolution of OLI in comparison with the common size of agricultural farms is its advantage for using in agricultural studies. OLI consist of nine spectral bands. The spatial resolution for bands 1 to 7 (coastal, blue, green, red, NIR, SWIR1, SWIR2) and band 9 (cirrus), is 30 m and for panchromatic band (band 8) is 15 m. A series of OLI images of the study area (path: 166, 167 & row: 35, 36) was acquired after wheat dormancy stage with no cloud or less than 10% amount of cloud in each year. The images were acquired in 2015 on April 4, April 20, April 29, May 6, May 31, June 7, June 16 and June 23 which there were coincided with tillering, stem extension (20th & 29th of April), booting, heading, dough development, physiological maturity and harvest ripe stages of wheat, respectively. Also, in 2016 the images on April 6, April 22, May 1, May 8, May 17, June 9, June 18 and June 25 were acquired which there were coincided with tillering, stem extension (22th of April & 1st of May), booting, heading, dough development, physiological maturity and harvest ripe stages of wheat, respectively.

The pre-processing of images was performed in ENVI 5.1 software. Since the OLI images had been corrected geometrically, we just did the second geo-referencing operation to ensure the correct location of terrain. The ground control points derived from 1:25000 topographic maps were applied for geometric correction. The geo-correction error less than 0.5 pixel (15 m) was attained. The nearest neighbor method and linear polynomial geometric model were used. Atmospheric corrections were performed using MODTRAN 4 model in FLAASH (Fast Line-of-sight Atmospheric Analysis of Spectral Hypercubes) module in the ENVI 5.1. The parameters used in the FLAASH were adjusted based on scene center location, sensor type, sensor altitude, ground elevation, pixel size, information about sensor flight date and weather conditions on the image acquisition date. The output of the FLAASH package was the surface reflectance of OLI images.

Extraction of vegetation indices

There are some studies which have suggested double or multi-band spectral indices to estimate bio-physical changes of crops (Garroutte et al., 2016; Li et al., 2016). This study included eight widely used VIs to investigate the strength of VIs during wheat gro­wing stages. These indices included: NDVI (normalized difference vegetation index), SAVI (soil adjusted vegetation index), EVI (enhanced vegetation index), EVI2 (2 bands enhanced vegetation index), NDWI (normalized difference water index), GreenNDVI (green normalized difference vegetation index), CVI (chlorophyll vegetation index) and CIgreen (green chlo­rophyll index). The VIs mentioned above can be classified in three groups based on their sensitivity to green biomass (NDVI, SAVI, EVI, EVI2), the liquid water content of vegetation (NDWI) and leaf chlorophyll (CVI, GreenNDVI, and CIgreen). Table 1 represents the formula of these indices.

Determining the best temporal stage for WAHD prediction

In each year, to determine the best temporal and phenological stage for predicting WAHD, correlation coefficients (R) between VIs and WAHD were calculated for different images acquired during the wheat growth phases. The date with the highest R values was recognized as the best temporal stage for WAHD prediction.

The artificial neural network model

The ANNs are flexible mathematical models that accomplish a computational simulation based on the behavior of human brain neurons. The ANN is a non-linear machine learning algorithm with a high potential for modeling and prediction (Foody, 2004). It is composed of artificial neuron groups which are interconnected with weighted links and could create a strong relation between inputs and outputs using a learning approach (Omer et al., 2016). The network architecture is the first step for developing an ANN model that is determined by artificial neurons and layers. A usual ANN consists of an input layer, an output layer, and several hidden layers. Also, each layer contains some neurons and activation functions. In this study, we used a Multi-layer perceptron (MLP) network which is extensively used type of ANNs in the community of remote sensing. The gradient descent with momentum (GDM) algorithm was used to train the network. For this algorithm, the weights are updated during the nth training iteration as follows (Omid et al., 2009):

Additionally, the set of weights are given by:

where W_ij⁽ⁿ⁾ is the weight between the jth neuron of the upper layer and the ith neuron of the lower layer, δj is the error of jth neuron, oi is the output value of the ith neuron of the previous layer, η is the learning rate and α is the momentum term. ΔW_ij⁽ⁿ⁾, is the gradient vector affiliated with the weights.

ANNs structure

To consider the efficacy of inputs (Green biomass sensitive, chlorophyll sensitive and water content sensitive VIs) on WAHD prediction, eight MLP type models were designed. Fig. 2 shows a typical scheme of an MLP model which comprises the input layer, the hidden layer, and the output layer. Different models were designed using various combination of VIs as the neurons in the input layer (ANN-1 to ANN-8). The input variables for ANN-1 to ANN-3 were green biomass sensitive VIs, chlorophyll sensitive VIs and water content sensitive VIs, respectively. The input variables for ANN-4 were EVI2 (green biomass sensitive), NDWI (water content sensitive) and GreenNDVI (chlorophyll sensitive) that each of them had the highest correlation with the observed WAHD in their groups. The input variables for ANN-5 was a combination of the NDVI and EVI2 (green biomass sensitive) and NDWI (water content sensitive). ANN-6 was a combination of green biomass sensitive VIs (NDVI, SAVI, EVI2) and chlorophyll sensitive VI (GreenNDVI). The input variables for ANN-7 were the combination of chlorophyll sensitive and water content sensitive VIs. The input variables for ANN-8 were the combination of green biomass sensitive and water content sensitive VIs (Table 2).

Figure 2. Scheme of the multi-layer perceptron (MLP) model structure used for wheat appropriate harvest date (WAHD).

For determining the optimum number of hidden layer neurons, the minimum value of RMSE was considered. In this study, the used activation function was hyperbolic tangent (TANH) for hidden and output layers. Also, the values of α=0.7 and η=0.1 were used.

The data set on 100 VIs for predicting WAHD was split into three parts: 70% of VIs for training, 15% for cross-validation and 15% for testing data. After sufficient training, the weights of the network were adapted and applied for validation to determine the model overall performance.

ANN model performance and validation

The performance of eight ANN constructed model using a different number of VIs was evaluated by statisti­cal parameters including coefficient of determination for prediction (test data) (R²_P), mean absolute error for prediction (MAEP), and the root mean square error for prediction (RMSEP). The values of R²_P, MAEP, and RMSEP were calculated using equations 3, 4 and 5, respectively.

where Ymea,i is the measured WAHD, Yest,i is the esti­mated WAHD, Ymea,i is the mean of measured values of WAHD, Yest,i is the mean of estimated values of WAHD, and N is the number of samples. The higher values for R²_P and lower values of MAEP and RMSEp re­present further precision and accuracy of the model. All calculations for ANN models were implemented using custom-written scripts of MATLAB R2015a software.

Harvest schedule regional (HSR) map

The predicted data of the best MLP model (in terms of precision and accuracy) was used for generating the HSR map using ArcGIS 10 software. Furthermore, the WAHD map was provided using the data of the test set. Moreover, both of the generated maps were compared with maps of measured data.

ResultsTop

The best temporal stage for WAHD prediction

Table 3 and Table 4 show the variation of R-values between VIs and observed WAHD in 2015 and 2016, respectively. The R-values were calculated to specify which temporal stage and VIs are appropriate for WAHD prediction. Also, for better comprehend the results in Tables 3 and 4, the temporal variations of R-values were presented in Fig. 3. Most of R-values had an ascendant trend before Jun 7, 2015 (Fig. 3a), and Jun 9, 2016 (Fig. 3b). However, a decline was observed in R-values after these dates. Based on the results the days of early of June were the best temporal phase for WAHD prediction which was coincided with dough development stage of wheat.

Figure 3. Temporal variations of R values between WAHD and VIs (a) in 2015 and (b) in 2016.

In 2015, the maximum R-values for NDVI, SAVI, EVI, EVI2, GreenNDVI, CIgreen, and NDWI were 0.729, 0.698, 0.744, 0.758, 0.594, 0.558 and 0.753 respectively that were obtained on Jun 7, but the R-value of CVI (0.563) was maximum on Jun 16 (Table 3).

As shown in Table 4, the maximum R-values were achieved on Jun 9 for NDVI (0.698), SAVI (0.686), EVI (0.709), EVI2 (0.763), and NDWI (0.747) in 2016. Additionally, the maximum of R-values for GreenNDVI, CVI, and CIgreen were 0.582, 0.554 and 0.541 on Jun 18. Considering the R-value as an indicator to assess the ability of VIs to estimate WAHD at wheat different growth stages, the June 7, 2015, and June 9, 2016, was selected for prediction.

Evaluation of ANN models

The results of designed models and the R²_P, MAEP, and RMSEP values for these models are summarized in Table 5. For these models, R²_P, MAEP and RMSEP values ranged from 0.611 to 0.925, 0.70 to 1.51, and 0.87 to 2.00, respectively. From the designed and trained networks, the ANN-5 presented better results than other networks. Therefore, the ANN-5 was selected as the best estimation model for WAHD (eq. 6):

where WAHDANN-5 is the predicted value of the WAHD using the ANN-5, and f(neti) is tanh (neti). Hence, the (1- 10-3) - MLP network, consists of three VIs as input variables (EVI2, NDWI, NDVI), ten neurons in the hidden layer and a single output variable (WAHD) was selected as the optimum network. The values of R²_P, MAEP, and RMSEP for this topology were obtained 0.925 (Fig. 4b), 0.70 and 0.87, respectively. Fig. 4a shows the comparison between predicted WAHD values (test data) with the measured data. It could be observed that the trends were similar and superimposed over the others in some parts. The second model that had good performance for WAHD estimation was ANN-8 (R²_P=0.882; MAEP=1.05; RMSEP=1.11). Three VIs (NDVI, EVI2, and NDWI) were common parameters for both models ANN-5 and ANN-8. In other words, NDVI, EVI2, and NDWI can predict WAHD better than other indices.

Figure 4. (a) Estimated (ANN-5) and measured WAHD on testing data. (b) Measured vs. estimated WAHD values (ANN5).

ANN-2 (R²_P=0.706; MAEP=1.51; RMSEP=2.00) which included GreenNDVI, CVI, and CIgreen had low accuracy for WAHD modeling. Therefore, it can be concluded that GreenNDVI, CVI, and CIgreen are not good indicators for WAHD prediction. In addition, ANN-3 model with individual input (NDWI) did not have high accuracy estimation for WAHD (R²_P=0.611; MAEP=1.41; RMSEP=1.86). But, it can be seen in ANN-7 model that adding NDWI to GreenNDVI, CVI, CIgreen VIs (ANN-2) led to better estimations for WAHD prediction (R²_P=0.759; MAEP=1.22; RMSEP=1.49). Besides, by adding the NDWI to vegetation indices of ANN-1, the accuracy of model improved ANN-8 (Table 5). In total, it is concluded that NDWI composition with other indicators could improve the forecast of WAHD.

HSR map

The HSR map using all data is shown in Fig. 5b, where ANN-5 model has been used for estimation as respects it had the best prediction accuracy (as explained in detail earlier). Also, the comparison between measured and estimated HSR maps is presented in Fig. 5, where the ANN-5 shows to have made good estimations of WAHD. Also, the trend of harvesting dates for different farms is the same in the two maps (Fig. 5).

Figure 5. Comparison between measured (a) and estimated (using ANN-5) (b) HSR maps for All data.

As shown in estimated and measured maps (Fig. 5), the WAHD of northern regions is after the southern regions. Therefore, the harvest operation should start from the farms in southeast regions in the maps and completed in northwest farms. The differences in sowing dates, topographic variations, and the angle of sunlight, temperature, precipitation and soil moisture amount are the reasons of WAHD variation among these farms.

The estimated (using test data of ANN-5) and measured WAHD maps are shown in Fig. 6a and 6b, respectively. As seen in Fig. 6, the ANN-5 could provide very good predictions of WAHD. However, there were few differences between measured and estimated WAHD, and one disagreement was labeled with a black circle in Fig. 6. Therefore, it can be concluded that the MLP model with 1-10-3 topology had very good performance for harvest date predictions.

Figure 6. Comparison between measured (a) and estimated (using ANN-5) (b) WAHD maps for Test data.

DiscussionTop

In this study, we present a method for predicting WAHD, based on Landsat-8 derived VIs and ANNs to optimize harvest campaign scheduling and farm management. The utilized approach is based on the spectral variation of wheat canopy during phenological stages. The VIs, using the combination of different spectral bands, could explain the vegetation conditions during wheat growth stages. The VIs used in this study were classified in three groups based on their sensitivity to green biomass (NDVI, SAVI, EVI, EVI2), the liquid water content of vegetation (NDWI) and leaf chlorophyll (CVI, GreenNDVI, and CIgreen).

In the present study, application of ANNs, as non-linear modeling techniques, had excellent performance for estimating WAHD. The performance and accuracy of MLP models to WAHD estimation was due to the key features of these type models: (i) the intrinsic abilities of ANNs, such as learning, cross-validating and flexible processing; (ii) generating non-linear patterns between inputs and outputs which lead to accurate estimations of complex and dynamic data; and (iii) the VIs included intensive data which correlate non-linearly with spatial based WAHD. Therefore, considering the nature of VIs data, the ANN modeling techniques are a very good substitute for linear statistical methods. As reported by Xie et al. (2009), the ANN models could provide an accurate estimation (R2=0.817 and RMSE=0.4236) for predicting above ground grassland biomass based on Landsat ETM+. Also, other researches which have been done using ANN techniques and VIs, suggest that utilizing of ANNs in remote sensing applications, has been very successful in the agriculture sector (Li et al., 2007; Fortin et al., 2010; Pantazi et al., 2016).

The ANN-5 had the best estimation of WAHD. The input variables of this model were EVI2, NDWI, and NDVI. EVI2 could be recognized as the best index for WAHD prediction because of its high R value with WAHD. EVI2 consists of NIR and red bands; strong contrast of leaf scattering in the NIR wavelength and chlorophyll absorption in the red wavelength, makes EVI2 very robust to explain the vegetation conditions. EVI2 not only has an improved sensitivity over high biomass, in comparison with the SAVI, but also minimizes soil influences (Jin et al., 2016). EVI2 can capture subtle changes in vegetation structure and condition, especially to discriminate between leaf area index surface and greenness for vegetation with various soil background reflectances (Rocha & Shaver, 2009). Therefore, EVI2 is a suitable indicator for WAHD prediction. Results achieved by Bolton & Friedl (2013) as well as by Wang et al. (2015), proved that EVI2 has a high ability to predict maize yield and estimate rice phenology, respectively. NDWI is sensitive to the total amounts of the liquid water content of vegetation canopies (Gao, 1996). The water content of wheat canopy varies during wheat growth stages, and NDWI could present these water content variations. As reported by Liu et al. (2006) and Bao et al. (2008), NDWI had a good correlation with wheat biophysical parameters and yield. Although in the present study NDWI did not have a good ability to WAHD estimation alone

(ANN-3 model; R²_P=0.611; MAEP=1.41; RMSEP=1.86), adding NDWI to other indices improved the estimation power of models (ANN-2 vs. ANN-7 and ANN1 vs. ANN8). NDVI is the most practical indicator which is widely used in vegetation monitoring. According to Suwannachatkul et al. (2014), the use of NDVI to estimate rice harvest date resulted in acceptable accuracy (about eight days). NDVI was the only indicator used for predicting. The results of the current study showed that the combination of VIs for WAHD prediction improved the estimation accuracy (about one day). This was also confirmed in the present study, where green biomass sensitive VIs (NDVI, SAVI, EVI, EVI2) could estimate WAHD better than chlorophyll sensitive VIs (CVI, GreenNDVI, CI green) by comparing models ANN-1 and ANN-2. Moreover, ANN-6 and ANN-7 represent this issue. Therefore, it is suggested that biomass-sensitive VIs be used in the harvesting date prediction researches.

The methodology applied in this paper is straightway and just needs satellite imagery bands to predict WAHD. The substantial advantage of the used method is its relative simplicity which makes it suitable for regional scale applications. The present method can be extended to other regions of the world. Also, other satellite sensors can be applied with this method. Moreover, the current model (ANN-5) may be helpful to regions with similar climate, where the sowing and harvesting dates are similar to the study region. The main limitation of this method is the accessibility of valid satellite imagery. The cloudy weather disrupts access to satellite images. A solution to this problem is the simultaneous use of several satellite images (Shang et al., 2014; Wang et al., 2015).

The current study presents a preliminary investigation of WAHD using satellite imagery. This method can be applied to other crops considering the crop conditions. Also, the combination of satellite remote sensing models, crop models, and meteorological statistical models will improve predictions.