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.

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 m ² . The size of plots chosen was 60×60 m ² 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.

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 (20 ^th & 29 ^th 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 (22 ^th of April & 1 ^st 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.

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 CI _green (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 CI _green ). Table 1 represents the formula of these indices.

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 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 n ^th training iteration as follows ( Omid et al. , 2009 ):

Additionally, the set of weights are given by:

where w_ij⁽ⁿ⁾ is the weight between the j ^th neuron of the upper layer and the i ^th neuron of the lower layer, δ _j is the error of j ^th neuron, o _i is the output value of the i ^th neuron of the previous layer, η is the learning rate and α is the momentum term. Δw_ij⁽ⁿ⁾, is the gradient vector affiliated with the weights.

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).

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.

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 (MAE _P ), and the root mean square error for prediction (RMSE _P ). The values of R_P², MAE _P , and RMSE _P were calculated using equations 3, 4 and 5, respectively.

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.

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.

In 2015, the maximum R -values for NDVI, SAVI, EVI, EVI2, GreenNDVI, CI _green , 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 CI _green 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.

The results of designed models and the R_P², MAE _P , and RMSE _P values for these models are summarized in Table 5. For these models, R_P², MAE _P and RMSE _P 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 WAHD _ANN-5 is the predicted value of the WAHD using the ANN-5, and f ( net _i ) is tanh ( net _i ). 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², MAE _P , and RMSE _P 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; MAE _P =1.05; RMSE _P =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.

ANN-2 (R_P²=0.706; MAE _P =1.51; RMSE _P =2.00) which included GreenNDVI, CVI, and CI _green had low accuracy for WAHD modeling. Therefore, it can be concluded that GreenNDVI, CVI, and CI _green 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; MAE _P =1.41; RMSE _P =1.86). But, it can be seen in ANN-7 model that adding NDWI to GreenNDVI, CVI, CI _green VIs (ANN-2) led to better estimations for WAHD prediction (R_P²=0.759; MAE _P =1.22; RMSE _P =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.

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).

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.