Monitoring of Wheat Growth Status and Mapping of Wheat Yield ’ s within-Field Spatial Variations Using Color Images Acquired from UAV-camera System

Abstract: Applications of remote sensing using unmanned aerial vehicle (UAV) in agriculture has proved to be an effective and efficient way of obtaining field information. In this study, we validated the feasibility of utilizing multi-temporal color images acquired from a low altitude UAV-camera system to monitor real-time wheat growth status and to map within-field spatial variations of wheat yield for smallholder wheat growers, which could serve as references for site-specific operations. Firstly, eight orthomosaic images covering a small winter wheat field were generated to monitor wheat growth status from heading stage to ripening stage in Hokkaido, Japan. Multi-temporal orthomosaic images indicated straightforward sense of canopy color changes and spatial variations of tiller densities. Besides, the last two orthomosaic images taken from about two weeks prior to harvesting also notified the occurrence of lodging by visual inspection, which could be used to generate navigation maps guiding drivers or autonomous harvesting vehicles to adjust operation speed according to specific lodging situations for less harvesting loss. Subsequently orthomosaic images were geo-referenced so that further study on stepwise regression analysis among nine wheat yield samples and five color vegetation indices (CVI) could be conducted, which showed that wheat yield correlated with four accumulative CVIs of visible-band difference vegetation index (VDVI), normalized green-blue difference index (NGBDI), green-red ratio index (GRRI), and excess green vegetation index (ExG), with the coefficient of determination and RMSE as 0.94 and 0.02, respectively. The average value of sampled wheat yield was 8.6 t/ha. The regression model was also validated by using leave-one-out cross validation (LOOCV) method, of which root-mean-square error of predication (RMSEP) was 0.06. Finally, based on the stepwise regression model, a map of estimated wheat yield was generated, so that within-field spatial variations of wheat yield, which was usually seen as general information on soil fertility, water potential, tiller density, etc., could be better understood for applications of site-specific or variable-rate operations. Average yield of the studied field was also calculated according to the map of wheat yield as 7.2 t/ha.


Introduction
Remote sensing has been successfully used as an effective method for obtaining field information through analysis of reflectance or radiance of specific bands' digital numbers [1,2].According to different types of platforms, agricultural remote sensing could generally be categorized into satellite remote sensing, aerial remote sensing, and near-ground remote sensing.Multi-spectral satellite Remote Sens. 2017, 9, 289 2 of 14 imagery has long been applied to detect vegetation areas, monitor crop growth status, and estimate crop yield, etc., in large scale [3,4].Aerial remote sensing using aircrafts or balloon has been introduced as supplementary method and often carried out as one-time operations.Aerial photography has been used to monitor crop growth status as regional or medium-scale applications of agricultural remote sensing ever since 1950s by using color or color-infrared cameras [5].Near-ground agricultural remote sensing is often referred to as frame-based or pillar-based applications [6].Recently, the cutting edge application of small fixed-and/or rotary-wing unmanned aerial vehicles (UAV) for spatial sampling or preliminary mapping for variable-rate operations has begun, since UAV imagery may be acquired more cost-effectively, with excellent maneuverability as well as increasing spatial resolution, and with greater safety when compared with manned aircraft [7][8][9].Agricultural application of UAV remote sensing by using color cameras instantly provides researchers and farmers with actual and intuitive visualization of crop growth status, since color images accentuate particular vegetation greenness and have been suggested to be less sensitive to variations of illumination conditions [10,11].Meanwhile, the application of color cameras also sharply decreases the high cost of remote sensing [12], since most digital cameras use a Bayer-pattern array of filters to obtain an RGB digital image, and the acquisition of near-infrared (NIR) band images usually requires an extra filter that converts digital numbers of either blue or red light in Bayer array into NIR readings through massive post-processing and calibration work [13].Rasmussen et al. investigated the reliability of four different vegetation indices derived from consumer-grade true color camera as well as a color-infrared camera that are mounted on UAVs for assessing experimental plots, and concluded that vegetation indices of UAV imagery have the same ability as ground-based recordings to quantify crop responses to experimental treatments, although such shortcomings like angular variations in reflectance, stitching, and ambient light fluctuation should be taken into consideration [14].Torres-Sánchez et al. mapped multi-temporal vegetation fraction in early-season wheat fields by using a UAV equipped with commercial color camera and studied the influence of flight altitude and days after sowing on the classification accuracy, which showed that visible spectral vegetation indices derived from low-altitude UAV-camera system could be used as a suitable tool to discriminate vegetation in wheat fields in the early season [15].Woebbecke et al. tested several color indices derived using chromatic coordinates and modified hue to distinguish vegetation from background such as bare soil and crop residues [16], among which the excess green vegetation index (ExG) provides a near-binary intensity image outlining a plant region of interest has been widely cited.Cui et al. evaluated the reliability of using color digital images to estimate above ground biomass at canopy level of winter wheat by taking pictures from one meter above the top of the wheat canopy [17].In short, most past agricultural studies on color images focused on individual level of crop or weed, and the point-source samplings of which are usually inevitably both time-consuming and have to be conducted under poor working condition.
Wheat (Triticum spp., or Triticum aestivum L.) has been among the most produced cereal grains for a long time.In 2009, global wheat production reached about 680 million tons, which made it the second most produced cereal grain after maize [18].In order to optimize wheat yield and grain quality, especially in terms of protein content that varies significantly depending on different agricultural practices, optimal agronomic management in accordance with wheat development stages is critical.It requires that farmers have a detailed understanding of wheat growth status during each specific development stage in wheat cultivation [19], which indicates that real-time monitoring of actual wheat growth status throughout the wheat growing season is of vital importance in helping wheat growers with management decision making.Meanwhile, it was frequently reported that reproductive growth of wheat after the flowering stage is closely related to grain yield, and many studies in recent years on wheat growth have indicated that accumulative NDVI values of multi-temporal satellite remote sensing images after flowering stage have good relationship with crop yield [2,20].Therefore, in this study, eight-color orthomosaic images acquired from a low altitude UAV-camera system were used to intuitively monitor and assess overall growth status of a winter wheat field from heading stage to ripening stage from early June to the end of July 2015 at the interval of about one week in Hokkaido, Japan.A multivariate analysis of field samples of wheat yield and accumulative CVIs derived from orthomosaic images was performed to estimate wheat yield.Subsequently, an original method of mapping wheat yields within-field spatial variations was proposed to provide reference for decision making in terms of site-specific agriculture based on the result of regression model.

Materials and Methods
The overall approach of this proposed method of monitoring wheat growth status as, well as mapping within-field wheat yield's spatial variations, consists of four key steps: acquisition of high resolution UAV photographs using the low altitude UAV-camera system; post-processing of UAV images including orthomosaicking, georeferencing, and radiometric normalization; extraction of color vegetation indices from post-processed othomosaic images; and mapping of with-field wheat yield's spatial variations through multivariate analysis between color vegetation indices and sampled grain yields, shown in Figure 1.This study adopted World Geodetic System 1984, or WGS84, as the coordinate system for geo-referenced images, maps, or any coordinates used in this paper if not particularly indicated.
Remote Sens. 2017, 9, 289 3 of 13 2015 at the interval of about one week in Hokkaido, Japan.A multivariate analysis of field samples of wheat yield and accumulative CVIs derived from orthomosaic images was performed to estimate wheat yield.Subsequently, an original method of mapping wheat yields within-field spatial variations was proposed to provide reference for decision making in terms of site-specific agriculture based on the result of regression model.

Materials and Methods
The overall approach of this proposed method of monitoring wheat growth status as, well as mapping within-field wheat yield's spatial variations, consists of four key steps: acquisition of high resolution UAV photographs using the low altitude UAV-camera system; post-processing of UAV images including orthomosaicking, georeferencing, and radiometric normalization; extraction of color vegetation indices from post-processed othomosaic images; and mapping of with-field wheat yield's spatial variations through multivariate analysis between color vegetation indices and sampled grain yields, shown in Figure 1.This study adopted World Geodetic System 1984, or WGS84, as the coordinate system for geo-referenced images, maps, or any coordinates used in this paper if not particularly indicated.

Figure 1.
Proposed method of estimating wheat yield by using unmanned aerial vehicle (UAV) images.

Field Site and Acquisition of UAV Images
Experiments were established on a winter wheat farmland located in Memuro, Hokkaido, Japan (around 42.902041°N-42.899607°Nand 142.977953°E-142.981734°E),shown in Figure 2. The lower left field (marked in black rectangle shown in Figure 3) was planted with the winter wheat variety of Kitahonami, occupying about 3.2 ha., which is the most widely planted winter wheat variety in Hokkaido and is reported to have taken up about 90% acreages of winter wheat in Hokkaido alone [21].The wheat field was planted around 25 September 2014 and harvested on 27 July 2015.Nine grain samples of Kitahonami were taken on 24 July 2015, three days prior to harvesting, and the position of each yield sample was measured by using Trimble SPS855 GNSS modular receiver in RTK-GPS mode with horizontal positioning accuracy of 8 mm [22].Regional annual precipitation of this area is around 957.3 mm, with average annual temperature of 6.1 °C [23].
In this experiment, a small quadrotor (ENROUTE CO., LTD., Fujimino, Japan), shown in Figure 4, was used as the platform of low altitude (about 100 m above ground level) UAV-camera system, and its performance parameters were listed in Table 1.A laptop installed with Ground Controlling Station (GCS) software was used to monitor and control the autonomous UAV flight through telemetry radio.Autonomous flights were conducted eight times at the interval of about one week from winter wheat's heading stage to ripening stage, on 2, 10, 19, and 25 June, 2015, and 2, 10, 16, and 24 July 2015 (at about 11:00 local time), using the flight paths that were designed beforehand in the GCS software, shown in Figure 4 as blue lines.
A SONY ILCE-6000 commercial digital camera was used to take pictures in continuous mode every two seconds (f/8, 1/500 sec, ISO 100), and the camera specification was also described in Table 1.During each flight, the camera was fixed on a two-axis gimbal, pointing vertically downwards

Field Site and Acquisition of UAV Images
Experiments were established on a winter wheat farmland located in Memuro, Hokkaido, Japan (around 42.902041 • N-42.899607 • N and 142.977953 • E-142.981734• E), shown in Figure 2. The lower left field (marked in black rectangle shown in Figure 3) was planted with the winter wheat variety of Kitahonami, occupying about 3.2 ha., which is the most widely planted winter wheat variety in Hokkaido and is reported to have taken up about 90% acreages of winter wheat in Hokkaido alone [21].The wheat field was planted around 25 September 2014 and harvested on 27 July 2015.Nine grain samples of Kitahonami were taken on 24 July 2015, three days prior to harvesting, and the position of each yield sample was measured by using Trimble SPS855 GNSS modular receiver in RTK-GPS mode with horizontal positioning accuracy of 8 mm [22].Regional annual precipitation of this area is around 957.3 mm, with average annual temperature of 6.1 • C [23].
In this experiment, a small quadrotor (ENROUTE CO., LTD., Fujimino, Japan), shown in Figure 4, was used as the platform of low altitude (about 100 m above ground level) UAV-camera system, and its performance parameters were listed in Table 1.A laptop installed with Ground Controlling Station (GCS) software was used to monitor and control the autonomous UAV flight through telemetry radio.Autonomous flights were conducted eight times at the interval of about one week from winter wheat's heading stage to ripening stage, on 2, 10, 19, and 25 June 2015, and 2, 10, 16, and 24 July 2015 (at about 11:00 local time), using the flight paths that were designed beforehand in the GCS software, shown in Figure 4 as blue lines.
and took about 120 photos covering two adjacent fields in order to get enough ground control points (GCPs) for georeferrencing the orthomosaic images in post-processing.The ground resolution of these pictures was about 2.5 cm.Every 120 individual photographs were stitched together as an orthomosaic image by using Agisoft Photoscan software (Agisoft LLC, Petersburg, Russia), which was shown in Figure 3.Each of the orthomosaic images have about 3600 × 2450 pixels in size, and the ground resolution reached up to about 12.5 cm after orthomosaicking.and took about 120 photos covering two adjacent fields in order to get enough ground control points (GCPs) for georeferrencing the orthomosaic images in post-processing.The ground resolution of these pictures was about 2.5 cm.Every 120 individual photographs were stitched together as an orthomosaic image by using Agisoft Photoscan software (Agisoft LLC, Petersburg, Russia), which was shown in Figure 3.Each of the orthomosaic images have about 3600 × 2450 pixels in size, and the ground resolution reached up to about 12.5 cm after orthomosaicking.and took about 120 photos covering two adjacent fields in order to get enough ground control points (GCPs) for georeferrencing the orthomosaic images in post-processing.The ground resolution of these pictures was about 2.5 cm.Every 120 individual photographs were stitched together as an orthomosaic image by using Agisoft Photoscan software (Agisoft LLC, Petersburg, Russia), which was shown in Figure 3.Each of the orthomosaic images have about 3600 × 2450 pixels in size, and the ground resolution reached up to about 12.5 cm after orthomosaicking.A SONY ILCE-6000 commercial digital camera was used to take pictures in continuous mode every two seconds (f/8, 1/500 sec, ISO 100), and the camera specification was also described in Table 1.During each flight, the camera was fixed on a two-axis gimbal, pointing vertically downwards and took about 120 photos covering two adjacent fields in order to get enough ground control points (GCPs) for georeferrencing the orthomosaic images in post-processing.The ground resolution of these pictures was about 2.5 cm.Every 120 individual photographs were stitched together as an orthomosaic image by using Agisoft Photoscan software (Agisoft LLC, Petersburg, Russia), which was shown in Figure 3.Each of the orthomosaic images have about 3600 × 2450 pixels in size, and the ground resolution reached up to about 12.5 cm after orthomosaicking.
Georeferencing, or geocoding, is the process of assigning geographic coordinates to data (usually an image file) that are spatial in nature but have no explicit geographic coordinates (for example aerial photographs).It is necessary to georeference such images so that thereafter they can be further studied using geographic information system (GIS) technology.The georeferencing process of the orthomosaic images was accomplished using ArcMap 10.2 (ESRI Inc., Redlands, AB, Canada) software by adopting the 1st order polynomial transformation method and taking eight wheat row corners that dispersed around the field as ground controlling points (GCP).The transformation created two least-square-fit equations by comparing the image space coordinates of the GCPs with the geographic coordinates (latitude and longitude) and translated the image coordinates of each pixel into geographic coordinates.These GCPs' geo-spatial coordinates were measured by using a Trimble SPS855 GNSS modular receiver in RTK-GPS mode.

Radiometric Normalization of Multi-Temporal UAV Images
Due to different illuminative situations, radiometric accuracy and consistency are difficult to maintain among multi-temporal remote sensing images.Relative radiometric correction, or radiometric normalization, was usually necessary to adjust multi-temporal remote sensing images to a set of reference data and compensate for the radiometric effects.In this paper, the pseudo-invariant features (PIF) method, which refers to ground objects whose reflectance values are nearly constant over time during a certain period [24], was used to perform radiometric normalization of these UAV orthomosaic images, band by band, respectively.Seven places along the road and five places on the roof were selected as PIFs in each othomosaic images.Around each PIF's location, mean value of pixels distributed within the area of about 0.25 m 2 were calculated in each orthomosaic image, so that the influence of abnormal values caused by foreign matters etc. could be decreased.Subsequently, according to the PIFs' values extracted from eight orthomosaic images taken on different dates, radiometric normalization models were built by performing linear regressions, taking the PIFs' values extracted from each orthomosaic image as predictive terms and reference data as response variable.The reference data were generated by averaging each PIF's multi-temporal pixel values of blue band, green band, and red band, respectively.The spatial distribution of the UAV images' PIFs was shown as dark dots on the roof and road in Figure 5.
images taken on different dates, radiometric normalization models were built by performing linear regressions, taking the PIFs' values extracted from each orthomosaic image as predictive terms and reference data as response variable.The reference data were generated by averaging each PIF's multi-temporal pixel values of blue band, green band, and red band, respectively.The spatial distribution of the UAV images' PIFs was shown as dark dots on the roof and road in Figure 5.

Field Sampling of Wheat Yield
Nine samples of wheat yield were measured by using a 1 m × 1 m square frame to separate samples of wheat canopies.Samples were selected after detailed visual inspection on eight orthomosaic images and field survey on the sampling day in order to representatively choose three samples of or nearby lodging area and six samples out of normal areas, of which the most flourished areas and sparse areas were both taken into consideration.Samples' spatial distribution was shown in Figure 6.The sampling operation was conducted at 24 July 2015, three days ahead of harvesting, by collecting wheat ears within the specified 1-square-meter section.These samples' geo-coordinates were acquired by using Trimble SPS855 GNSS modular receiver in RTK-GPS mode.After threshing, grain weight of each sample was calculated and converted to 12.5% moisture, listed in Table 2.
Remote Sens. 2017, 9, 289 6 of 13 Nine samples of wheat yield were measured by using a 1 m × 1 m square frame to separate samples of wheat canopies.Samples were selected after detailed visual inspection on eight orthomosaic images and field survey on the sampling day in order to representatively choose three samples of or nearby lodging area and six samples out of normal areas, of which the most flourished areas and sparse areas were both taken into consideration.Samples' spatial distribution was shown in Figure 6.The sampling operation was conducted at 24 July 2015, three days ahead of harvesting, by collecting wheat ears within the specified 1-square-meter section.These samples' geo-coordinates were acquired by using Trimble SPS855 GNSS modular receiver in RTK-GPS mode.After threshing, grain weight of each sample was calculated and converted to 12.5% moisture, listed in Table 2.

Generating CVI Maps Based on UAV Images
Vegetation index map refers to as a scalar image in which each pixel has only one single brightness value.The pixel values are often calculated from reflectance or radiance of specific bands of remote sensing images.In recent years, several CVIs based on color images that are different from the NDVI associated with near-infrared band were proposed to identify vegetative features such as ExG mentioned above.Other CVIs were also introduced in this study, including visible-band difference vegetation index (VDVI) [11], normalized green-red difference index (NGRDI) [25], normalized green-blue difference index (NGBDI) [11], green-red ratio index (GRRI)

Generating CVI Maps Based on UAV Images
Vegetation index map refers to as a scalar image in which each pixel has only one single brightness value.The pixel values are often calculated from reflectance or radiance of specific bands of remote sensing images.In recent years, several CVIs based on color images that are different from the NDVI associated with near-infrared band were proposed to identify vegetative features such as ExG mentioned above.Other CVIs were also introduced in this study, including visible-band difference vegetation index (VDVI) [11], normalized green-red difference index (NGRDI) [25], normalized green-blue difference index (NGBDI) [11], green-red ratio index (GRRI) [26], and ExG [16], which were expressed respectively in the following equations: where B, G, R denotes the radiometric normalized pixel values of each orthomosaic images' blue, green, and red band, respectively.Accordingly, orthomosaic images taken on eight different dates were all used to generate CVI maps of ExG, NGBDI, GRRI, NGRDI, and VDVI as scalar images by using ENVI software (Exelis VIS, Inc., Boulder, CO, USA), in which each pixel was determined from its existing multi-band radiometric normalized pixel values and using band math functions.Figure 7 illustrates examples of CVI maps based on orthomosaic images taken on 2 June 2015, as well as the accumulative ExG map, whilst CVI maps based on other orthomosaic images and accumulative CVI maps were omitted here due to limited space.The accumulative CVI maps were acquired, respectively, by overlapping each corresponding CVI maps and adding pixel values of these scalar images on different dates point for point.After applying a mean filter of 7 × 7 pixels (which covers about 1-square-meter area) to the accumulative CVI maps, brightness values of pixels that have the same geo-spatial coordinates with the nine wheat samples were extracted out of each accumulative CVI map.The extracted accumulative CVI values were used to conduct stepwise regression analysis with sampled wheat yield data.Then, according to the regression model among wheat yield data and values of accumulative CVIs, yield map was generated accordingly.different dates point for point.After applying a mean filter of 7 × 7 pixels (which covers about 1-square-meter area) to the accumulative CVI maps, brightness values of pixels that have the same geo-spatial coordinates with the nine wheat samples were extracted out of each accumulative CVI map.The extracted accumulative CVI values were used to conduct stepwise regression analysis with sampled wheat yield data.Then, according to the regression model among wheat yield data and values of accumulative CVIs, yield map was generated accordingly.

Monitoring of Wheat Growth Status
From multi-temporal UAV orthomosaic images (shown in Figure 8), we can get a straightforward visualization of the rapid change of wheat growth status through image

Monitoring of Wheat Growth Status
From multi-temporal UAV orthomosaic images (shown in Figure 8), we can get a straightforward visualization of the rapid change of wheat growth status through image interpretation.We can see that canopy greenness reached peak condition on 10 July 2015 and the process of yellowing began thereafter.We can also see within-field variation of wheat tiller densities, especially from early stage of wheat growth, from the image taken on 2 June 2015.To be specific, the areas circled in red had relatively higher level of tiller densities, which compounded with other environmental influences such as rainfall, wind, etc. and caused occurrence of lodging before harvesting when over-luxuriant canopies failed to support heavy wheat ears (see images from 16 July and 24 July, where these areas are circled in black).Figure 9 showed a close-shot photograph of the lodging spot in test wheat field, taken on 16 July 2015.The orthomosaic UAV images taken at early stage of wheat growth could be used as prescription of variable-rate fertilizing to avoid, or alleviate the occurrence of lodging by reducing dozes of fertilizer around over-high tiller density areas; while the orthomosaic image taken ahead of harvesting could be practically served as references that guide either drivers of combine harvesters or autonomous harvesting vehicles to adjust operation speed according to the lodging situations, since lodging has been widely considered as the main cause of the deteriorating of grain quality and high loss rate of harvesting.early stage of wheat growth could be used as prescription of variable-rate fertilizing to avoid, or alleviate the occurrence of lodging by reducing dozes of fertilizer around over-high tiller density areas; while the orthomosaic image taken ahead of harvesting could be practically served as references that guide either drivers of combine harvesters or autonomous harvesting vehicles to adjust operation speed according to the lodging situations, since lodging has been widely considered as the main cause of the deteriorating of grain quality and high loss rate of harvesting.

Evaluation of Radiometric Normalization of Multi-Temporal Orthomosaic Images
According to PIFs' averaged band values, as well as reference data, radiometric normalization models of orthomosaic images taken on different dates, band by band, were built and shown in Table 3. From the R-squared values of regression models, we can see that the correlation between PIFs' blue band and reference data showed the most significant irrelevance, while the green band and red band expressed consistently higher relevance and less variation, indicating that the blue band was more susceptible to influences of different photographing conditions.We can also conclude that the image of 10 June 2015 was taken under the most deviated photographing condition when compared to other images.Based on linear-regression normalization models, each band of all images was normalized to the reference data and the effects caused by different photographing condition could be compensated.

Evaluation of Radiometric Normalization of Multi-Temporal Orthomosaic Images
According to PIFs' averaged band values, as well as reference data, radiometric normalization models of orthomosaic images taken on different dates, band by band, were built and shown in Table 3. From the R-squared values of regression models, we can see that the correlation between PIFs' blue band and reference data showed the most significant irrelevance, while the green band and red band expressed consistently higher relevance and less variation, indicating that the blue band was more susceptible to influences of different photographing conditions.We can also conclude that the image of 10 June 2015 was taken under the most deviated photographing condition when compared to other images.Based on linear-regression normalization models, each band of all images was normalized to the reference data and the effects caused by different photographing condition could be compensated.

Mapping of Wheat Yield's within-Field Spatial Variations
As mentioned in Section 1, it was reported that accumulative CVI values of multi-temporal remote sensing images after flowering stage have good relationship with crop yield.In this paper, we conducted stepwise regression analysis of sampled wheat yield with five different accumulative CVIs extracted from eight orthomosaic images that covered a winter wheat field from heading stage to ripening stage.As described in Sections 2 and 4, accumulative values of each CVI were obtained by extracting averaged values of pixels that fall within 1-square-meter area around the locations of nine wheat yield samples in each accumulative CVI maps of ExG, NGBDI, GRRI, NGRDI, and VDVI, listed in Table 4.By using stepwise method, a regression analysis was performed in MATLAB R2013a (The MathWorks, Inc., Natick, MA, USA) among the response variable of sampled wheat yield (listed in Table 2) and the predictive variables of accumulative CVI s (listed in Table 4).The result showed that the variable of NGRDI was removed from the stepwise regression model, while the rest variable of CVIs were included to fit the regression model expressed as following, with coefficient of determination as 0.94 and RMSE = 0.02: Y = −6.19− 6.78 × X1 + 3.45 × X3 + 0.88 × X4 + 0.003 × X5 (6) where Y, X1, X3, X4, and X5 denotes estimated wheat yield, accumulative VDVI, NGBDI, GRRI, and ExG, respectively.Subsequently, leave-one-out cross-validation (LOOCV) [27] was also conducted in MATLAB by building nine linear regression models which uses eight set variables of sampled wheat yield and values of CVIs of VDVI, NGBDI, GRRI, and ExG as training data, whilst leaves one set of variables as test data.According to Equation (7) [28], root-mean-square error of prediction (RMSEP) and correlation coefficient r were calculated as 0.06 and 0.69, respectively.
where y p and y m denotes predicted value of wheat yield according to each linear regression models by using training data mentioned above and test data of measured grain weight per square meter, respectively.Based on the regression model expressed as Equation ( 6), wheat yield was calculated by extracting each pixel values of accumulative CVI maps of VDVI, NGBDI, GRRI, and ExG in ENVI software, and map of wheat yield was generated accordingly.Wheat yield's within-field spatial variations could be observed from the map shown in Figure 10.From wheat yield map we obtained the information that about 25.8% areas in the studied field had grain weight per square meter below 0.5 kg; whilst grain weight per square meter of most areas reached between 0.5-1.5 kg, occupying about 50.4% acreages; and the mean value of grain weight per square meter was calculated as 0.72 kg, which indicated the estimated average yield of the studied field as 7.2 t/ha.

Uncertainties, Errors, and Accuracies
The civilian use of UAVs opens up a new way of obtaining field information that is both time and cost efficiently for precision agriculture.Wang [11] et al. successfully differentiated vegetation areas from non-vegetation areas by analyzing color images acquired from a UAV.Rasmussen [14] et al. investigated four different vegetation indices acquired from a color camera and a color-infrared camera by using both a fixed-wing UAV and a rotary-wing UAV, and concluded that vegetation indices based on UAV imagery have the same ability to quantify crop responses with ground-based recordings.However, when UAV imagery was applied into quantitative remote sensing, special attention should be paid to the process of ground truth samplings.In this study, we carefully selected nine grain weight samples to conduct stepwise regression analysis with five color vegetation indices, and cross validation also showed good predicative validity.The small size of the wheat field might also contribute to a certain extent to the significant validity of the regression model, and we cannot guarantee the applicability into other crop field before massive sampling is conducted over several farmlands around a large area.Since the altitude of UAV flight affects image resolution severely, which in turn changes the vegetation index's value by weakening or strengthening background (soil or crop residues) interferences, appropriate flight altitude should also be taken into consideration when conducting quantitative analysis.Study on vegetation indices also showed that accumulative VDVI+NGBDI+ExG correlated with accumulative NGRDI+GRRI with the correlation coefficient as −0.84, whist accumulative VDVI correlated with sampled grain weight per square meter with the correlation coefficient as 0.85.Therefore, regression analysis using different combination of vegetation indices might also affect the validity of regression model.

Conclusions
From multi-temporal UAV orthomosaic images we could visualize the rapid change of wheat growth status through image interpretation and discerned that the canopy greenness of the wheat field of study reached peak condition on 2 July 2015 and the process of yellowing began thereafter.We could also see within-field variation of wheat tiller densities, especially from an early stage of wheat growth from the image taken on 2 June 2015.The occurrence of lodging could also be spotted in orthomosaic images taken on 16 July and 24 July 2015.Therefore, we reached the conclusion that orthomosaic UAV images taken at early stage of wheat growth could be used as prescription of variable-rate fertilizing to avoid, or alleviate, the occurrence of lodging by reducing dozes of fertilizer around over-high tiller density areas.The orthomosaic image taken ahead of harvesting could be used to generate navigation maps that guide either drivers of combine harvesters, or autonomous harvesting vehicles, to adjust operation speed according to the specific lodging situations.

Uncertainties, Errors, and Accuracies
The civilian use of UAVs opens up a new way of obtaining field information that is both time and cost efficiently for precision agriculture.Wang et al. [11] successfully differentiated vegetation areas from non-vegetation areas by analyzing color images acquired from a UAV.Rasmussen et al. [14] investigated four different vegetation indices acquired from a color camera and a color-infrared camera by using both a fixed-wing UAV and a rotary-wing UAV, and concluded that vegetation indices based on UAV imagery have the same ability to quantify crop responses with ground-based recordings.However, when UAV imagery was applied into quantitative remote sensing, special attention should be paid to the process of ground truth samplings.In this study, we carefully selected nine grain weight samples to conduct stepwise regression analysis with five color vegetation indices, and cross validation also showed good predicative validity.The small size of the wheat field might also contribute to a certain extent to the significant validity of the regression model, and we cannot guarantee the applicability into other crop field before massive sampling is conducted over several farmlands around a large area.Since the altitude of UAV flight affects image resolution severely, which in turn changes the vegetation index's value by weakening or strengthening background (soil or crop residues) interferences, appropriate flight altitude should also be taken into consideration when conducting quantitative analysis.Study on vegetation indices also showed that accumulative VDVI + NGBDI + ExG correlated with accumulative NGRDI+GRRI with the correlation coefficient as −0.84, whist accumulative VDVI correlated with sampled grain weight per square meter with the correlation coefficient as 0.85.Therefore, regression analysis using different combination of vegetation indices might also affect the validity of regression model.

Conclusions
From multi-temporal UAV orthomosaic images we could visualize the rapid change of wheat growth status through image interpretation and discerned that the canopy greenness of the wheat field of study reached peak condition on 2 July 2015 and the process of yellowing began thereafter.We could also see within-field variation of wheat tiller densities, especially from an early stage of wheat growth from the image taken on 2 June 2015.The occurrence of lodging could also be spotted in orthomosaic images taken on 16 July and 24 July 2015.Therefore, we reached the conclusion that orthomosaic UAV images taken at early stage of wheat growth could be used as prescription of variable-rate fertilizing to avoid, or alleviate, the occurrence of lodging by reducing dozes of fertilizer around over-high tiller density areas.The orthomosaic image taken ahead of harvesting could be used to generate navigation maps that guide either drivers of combine harvesters, or autonomous harvesting vehicles, to adjust operation speed according to the specific lodging situations.
Through stepwise regression analysis of the response variable of sampled grain weight per square meter and the predictive variables of accumulative color vegetation indices, we can conclude that only the variable of normalized green-red difference index was removed from the stepwise regression model due to insignificant p-value, whilst the rest variables of visible-band difference vegetation index (the normalized green-blue difference index, green-red ratio index, and excess green vegetation index) were included to fit the regression model, with coefficient of determination and RMSE as 0.94 and 0.02, respectively.The averaged value of sampled grain weight per square meter was 0.86 kg.The regression model was also validated by using leave-one-out cross validation method, which showed that the root-mean-square error of predication of the regression model was 0.06.
Based on the stepwise regression model, a map of estimated grain weight per square meter (yield map) was generated by extracting each pixel values out of the maps of accumulative vegetation indices.Within-field spatial variations of wheat yield could be observed from the map, which could be seen as the comprehensive presentation of the spatial variations of soil fertility, tiller density, effective water potential, canopy aeration condition, and so on.We also obtained general information of the studied field that about 25.8% areas in the studied field had grain weight per square meter below 0.5 kg, whilst grain weight per square meter of most areas reached between 0.5-1.5 kg, occupying about 50.4% acreages.The mean value of grain weight per square meter was calculated as 0.72 kg, which indicated the estimated average yield of the studied field as 7.2 t/ha.

Figure 1 .
Figure 1.Proposed method of estimating wheat yield by using unmanned aerial vehicle (UAV) images.

Figure 2 .
Figure 2. Field site of the test farmland in Memuro, Hokkaido, Japan.

Figure 3 .
Figure 3. Winter wheat field of study was shown in black rectangle (ground control points were marked as black dots).

Figure 2 .
Figure 2. Field site of the test farmland in Memuro, Hokkaido, Japan.

Figure 3 .
Figure 3. Winter wheat field of study was shown in black rectangle (ground control points were marked as black dots).

Figure 2 .
Figure 2. Field site of the test farmland in Memuro, Hokkaido, Japan.

Figure 3 .
Figure 3. Winter wheat field of study was shown in black rectangle (ground control points were marked as black dots).

Figure 7 .
Figure 7. CVI maps on 2 June 2015 and accumulative ExG map of the test wheat field.

Figure 8 .
Figure 8. Orthomosaic images of the test wheat field from heading stage to harvesting.

142°58Figure 8 .
Figure 8. Orthomosaic images of the test wheat field from heading stage to harvesting.

Figure 9 .
Figure 9. Close-shot photograph of the lodging spot in test wheat field, taken on 16 July 2015.

Figure 9 .
Figure 9. Close-shot photograph of the lodging spot in test wheat field, taken on 16 July 2015.

Figure 10 .
Figure 10.Map of wheat yield (expressed as grain weight per square meter).

Figure 10 .
Figure 10.Map of wheat yield (expressed as grain weight per square meter).

Table 1 .
Performance parameters of UAV-camera system.

Table 2 .
Samples of wheat yield.

Table 2 .
Samples of wheat yield.

Table 4 .
Accumulative values of each color vegetation index.