Simulation of Reﬂectance and Vegetation Indices for Unmanned Aerial Vehicle (UAV) Monitoring of Paddy Fields

: Reﬂectance and vegetation indices obtained from aerial images are often used for monitoring crop ﬁelds. In recent years, unmanned aerial vehicles (UAVs) have become popular and aerial images have been collected under various solar radiation conditions. The value of observed reﬂectance and vegetation indices are considered to be a ﬀ ected by solar radiation conditions, which may lead to inaccurate estimations of crop growth. In this study, in order to evaluate the e ﬀ ect of solar radiation conditions on aerial images, canopy reﬂectance in paddy ﬁelds was simulated by a radiative transfer model, FLiES (Forest Light Environmental Simulator), for various solar radiation conditions and canopy structures. Several parameters including solar zenith angle, proportion of di ﬀ use light for incident sunlight, plant height, coordinates of plants and leaf area density, were tested in FLiES. The simulation results showed that the solar zenith angle did not vary the canopy reﬂectance under the conditions of the proportion of di ﬀ use light at 1.0, but the variation was greater at lower proportions of di ﬀ use light. The di ﬀ erence in reﬂectance caused by solar radiation was 0.01 and 0.1 at the maximum for red and near-infrared bands, respectively. The simulation results also showed that the di ﬀ erences in reﬂectance a ﬀ ect vegetation indices (Simple Ratio (SR), Normalized Di ﬀ erence Vegetation Index (NDVI) and Enhanced Vegetation Index 2 (EVI2)). The variation caused by solar radiation conditions was the least for NDVI and the greatest for SR. However, NDVI was saturated at the least leaf area index (LAI), whereas SR was only slightly saturated. EVI2 was intermediate between SR and NDVI, both in terms of variation and saturation. The simulated reﬂectance and vegetation indices were similar to those obtained from the aerial images collected in the farmers’ paddy ﬁelds. These results suggest that a large proportion of di ﬀ use light (close to 1.0) or a middle range of solar zenith angle (45 to 65 degrees) may be desirable for UAV monitoring. However, to maintain ﬂexibility of time and occasion for UAV monitoring, EVI2 should be used to evaluate crop growth, although calibration based on solar radiation conditions is recommended.


Introduction
An appropriate cultivation management strategy in crop fields should be employed to produce crop grain with high quality and quantity. The appropriate cultivation management first requires monitoring crop growth for farmers to decide suitable management, such as fertilization and water control. The leaf area index (LAI) is often used as an index for monitoring crop growth in crop fields because LAI is an appearance trait of crop status and is important for crop production [1][2][3][4]. LAI could be measured manually at the ground level using measuring instruments, such as a plant canopy analyzer [5]. The LAI measurement in crop fields is often restricted by labor and time because the measurement can only cover a limited range.
For this restriction, several efforts to estimate LAI using remote sensing techniques have been examined in past studies [6][7][8]. The estimation methods in these studies mostly used vegetation indices calculated using the canopy reflectance obtained from multispectral images taken in crop fields. For this purpose, various types of vegetation indices have been developed, such as Normalized Difference Vegetation Index (NDVI) and Soil-Adjusted Vegetation Index (SAVI) [9][10][11][12][13][14]. For example, Ali et al. estimated wheat LAI based on NDVI and SAVI and improved the vegetation indices by adding a red-edge spectral band [6]. These vegetation indices are often calculated from satellite remote sensing data, however, end users, such as farmers, have difficulties using the data because of the cost, especially with high spatial resolution, and the inflexibility of employing remote sensing techniques.
In recent years, since access to aerial images is becoming relatively easier as unmanned aerial vehicles (UAVs, such as drones) become more common [15], the utilization of aerial images by farmers in their crop fields has been expected. Some efforts to apply remote sensing techniques to aerial images taken by UAVs have also been conducted. Information obtained by analyzing aerial images taken by UAV was utilized for breeding [16,17], cultivation management [18,19] and yield estimation [20,21]. For example, Mukoyama et al. attempted to estimate the chlorophyll gauge value of rice by applying the vegetation index using the green and near-infrared spectral bands to the hyperspectral image [19]. Maki et al. evaluated the growth of LAI from multitemporal aerial images to estimate rice yield in agricultural fields [21].
UAV can fly at any time, except during restricted conditions, such as bad weather. Therefore, compared with satellite images, a UAV image can be obtained under various solar radiation conditions including various proportions of diffuse light for incident sunlight and various incident angles of direct sunlight (solar zenith angle). Different solar radiation conditions cause a difference in the observed reflectance for the same object [22,23]. Consequently, the difference in observed reflectance might affect the vegetation index, even though the vegetation index adjusts for such difference by using ratio or normalization of reflectance values. The vegetation index influenced by solar radiation conditions could produce inaccurate crop growth status and mislead farmers. However, few studies take the effect of solar radiation conditions into consideration when using a vegetation index for UAV monitoring of crop fields.
This study simulated the reflectance of a crop canopy with various solar radiation conditions using a radiative transfer model to evaluate the effect of solar radiation conditions on canopy reflectance and vegetation indices. Paddy fields were selected because rice is a major staple food and occupies a large portion of cultivated fields in Asia. To evaluate the effect of solar radiation conditions at various growth periods, the simulation was conducted by changing the parameters of canopy structure. For example, LAI and plant height corresponded with rice growth. The simulated reflectance and vegetation indices were plotted against LAI to evaluate the variation against solar radiation conditions. To validate the simulation results, aerial multispectral images were taken by UAV during several growth periods in farmers' paddy fields, and the canopy reflectance obtained from the images was compared with the simulated reflectance. Based on these results, optimal solar radiation conditions and vegetation indices for UAV monitoring of crop fields were discussed.

Simulation of Canopy Reflectance by a Radiative Transfer Model
A radiative transfer model was used as a method to simulate canopy reflectance. Photons that are transferred from the sun to the atmosphere to reach the observer, such as a spectral sensor, were simulated through scattering by atmospheric particles and reflection, transmission and absorption by the vegetation canopy. These processes are modeled by a radiative transfer model. Canopy reflectance is calculated by the ratio of the photons at the observer to those inputted to the canopy. In this study, one of the representative radiative transfer models, FLiES (Forest Light Environmental Simulator [22]), was applied to calculate canopy reflectance. FLiES simulates radiative transfer in the atmosphere and canopy based on the Monte Carlo ray tracing technique and calculates various plant light environmental conditions, including the bidirectional reflectance factor, absorbed photosynthetically active radiation and irradiance at the top and bottom of the canopy. FLiES was originally developed to simulate canopy reflectance in forests and has been applied in past studies [24][25][26][27][28][29], but it has also been used to simulate annual crops in previous studies [14,30].
FLiES allows users to evaluate the effects of parameters characterizing solar radiation conditions and canopy structure on canopy reflectance. The parameters used in this study included the proportion of diffuse light for incident sunlight and solar zenith angle for solar radiation conditions, plant coordinates, plant height, leaf reflectance, soil reflectance and leaf area density (LAD) for canopy structure (Table 1). Canopy LAI is calculated by multiplying the number of plants per unit area, the plant foliage volume and the LAD of the foliage. In this study, ellipsoid was applied as the foliage shape, and its volume was calculated using plant height and diameter that was assumed to correspond to hill interval. Canopy reflectance was simulated for two spectral bands: red (640-680 nm; RED) and near-infrared (770-810 nm; NIR). Three vegetation indices defined using these two bands were validated in this study: Simple Ratio (SR) [11], Enhanced Vegetation Index 2 (EVI2) [13] and Normalized Difference Vegetation Index (NDVI) [9]. These vegetation indices are defined as Equations (1), (2) and (3), respectively.

Verification of Simulated Reflectance by the Field Measurement
In this study, to verify the canopy reflectance by FLiES, aerial multispectral images were collected by UAV at several growth periods (20 June, 6 July, and 2 August, 2017) in farmers' paddy fields in Sendai, Miyagi prefecture, Japan (140 • 58'14"E, 38 • 13'2"N, Figure 1). Four farmers' paddy fields where cultivar 'Hitomebore', the most popular cultivar in Miyagi prefecture, was planted were selected. The acreages of the field were 0.84 (Field 1), 0.86 (Field 2), 0.9 (Field 3), 0.88 ha (Field 4), and the total acreage was 3.48 ha. The distance between planting rows was 0.3 m, and the hill interval in a row was 0.2 m. Twenty points in each field, for a total of eighty points, were used for field measurement and verification. Location and the orthomosaic reflectance image made from the aerial multispectral images taken by unmanned aerial vehicle (UAV) in the measured paddy fields. The imagery on the left side was taken from the Geospatial Information Authority of Japan [31]. The imagery on the right side was the orthomosaic reflectance image (6 July 2017), which was the false color RGB composite image, allocated a green spectral band to B, a red (RED) spectral band to G, and a near-infrared (NIR) spectral band to R.
The multispectral imaging sensor 'Sequoia' (Parrot, Paris, France) mounted on the UAV 'Solo' (3D Robotics, Berkeley, USA) was used to collect the aerial multispectral images. An aerial multispectral image can cover 0.4 ha at 60 m of altitude, and approximately 400 images were taken to cover the fields. All aerial multispectral images were taken at a 70% overlap rate both for the forward lap and side laps, and adjusted to create the orthomosaic reflectance image, with reflectance serving as a pixel value, using the image processing software 'Pix4D mapper' (Pix4D, Lausanne, Switzerland).
The Pix4D mapper extracted more than 10,000 key points, which were specific features in the multispectral image and were used to match with the other multispectral images. Based on the extracted key points, automatic aerial triangulation and bundle block adjustment were conducted to create a point cloud. The mean reprojection errors of the block bundle adjustment, which should be less than one pixel [32], were less than 0.3 pixels. A camera model optimization was also conducted to calibrate the internal and external parameters of the imaging sensor. The relative difference between the initial and the optimized parameters, which indicates the quality of camera model optimization and should be less than 5% [32], was less than 1%. The orthomosaic reflectance image with a spatial resolution of 0.05 m was created using the digital surface model (DSM). The reflectance was corrected by a radiometric calibration with the images of the calibration target (AIRINOV, Paris, France). These images were taken at the ground level at both the start and end of the image collection by UAV. The radiometric calibration also used information including sunlight irradiance, ISO, aperture, and shutter speed, which were recorded by imaging sensors and sunshine sensors mounted on UAV. Thirty-two ground control points were set at the corner of each field and uniformly distributed in the image. The points were used to georeference the orthomosaic reflectance images, and those coordinates were obtained from the aerial ortho image published as the digital Japan basic map by the Geospatial Information Authority of Japan (GSI) [31]. The aerial ortho image of the GSI was georeferenced at a location accuracy of less than 1 m [33]. Considering this location accuracy of the aerial ortho image and the locational deviation of the manual measurement of LAI, the average reflectance in a 3 m-diameter circle, centered on the field measurement point, was obtained from the orthomosaic reflectance image to compare with the canopy reflectance simulated by FLiES.
During the same growth period that the images were collected, plant height and LAD were measured at each of the eighty points. Plant height was measured by a scale. LAD was calculated by dividing LAI by the number of plants and the foliage volume. The foliage volume was calculated as the ellipsoid volume using the measured plant height and the diameter, which was assumed to correspond to the hill interval in paddy fields. LAI was measured using the plant canopy analyzer 'LAI-2200' (LI-COR, Lincoln, USA) [5]. The solar radiation conditions, proportion of diffuse light and solar zenith angle at the aerial multispectral images collected by UAV were 0.23 and 54.4 degrees on 20 June, 0.34 and 23.6 degrees on 6 July, and 0.98 and 64.6 degrees on 2 August, respectively. These solar radiation conditions were assumed to be constant during the image collection, because each UAV flight was conducted under relatively stable conditions with respect to the solar radiation condition. The proportion of diffuse light for incident sunlight was estimated by the method proposed by Erbs et al. [34]. The solar radiation measured by the Japan Meteorological Agency [35] was used in the method. The solar zenith angle was calculated based on the coordinates of the measured paddy fields and image collection date. Based on the measured plant height, LAD, and the solar radiation conditions in the paddy fields, the canopy reflectance was simulated by FLiES for each measurement point. Figure 2 shows the results of the canopy reflectance simulation for a paddy field using FLiES for the combination of solar radiation conditions and canopy structure parameters shown in Table 1. The simulated reflectance values were plotted against LAI, which was applied as the target parameter for UAV monitoring in this study.

Simulation of Reflectance Values
Although the solar zenith angle did not affect the reflectance of either RED or NIR with the proportion of diffuse light at 1.0, the effect was greater with a lower proportion of diffuse light. In the case of the proportion of diffuse light at 0.0 and 0.5, the small solar zenith angle (5 and 25 degrees) increased the reflectance of RED, whereas the large zenith angle (85 degrees) increased the reflectance of NIR. The greatest difference between the minimum and the maximum reflectance caused by solar radiation conditions for NIR was 0.1 at 0.5 of LAI, and that for RED was 0.01 at 4.   Figure 3 shows the results of the vegetation indices calculated with the reflectance in Figure 2. The vegetation indices with solar zenith angle of 5 and 85 degrees were excluded from the figure. The vegetation indices were plotted against LAI, which was applied as the target parameter for UAV monitoring in this study. The vegetation indices were the highest with conditions of the proportion of diffuse light at 1.0 and lowest with conditions of the proportion of diffuse light at 0.0. The solar zenith angle at 25 degrees produced relatively lower vegetation indices than that at 45 and 65 degrees, but the effect of solar zenith angle on canopy reflectance was relatively less than that of the proportion of diffuse light. The difference between the greatest and the least vegetation indices increased as LAI increased. The difference for SR was the greatest, and that for NDVI was the least. However, NDVI showed a saturation of index value at lower LAI than the other indices. Figure 4 shows the comparison of the average of reflectance obtained from the reflectance image made from the aerial multispectral images and the reflectance simulated by FLiES based on the parameters obtained in the field measurement.

Comparison of the Field Measurement Data and the Simulation Results
Although the standard deviations of reflectance obtained from the reflectance image were greater than the simulated reflectance, their averages were similar ( Table 2). For RED, the simulated reflectance was slightly underestimated on 20 June and slightly overestimated on 6 July and 2 August for the reflectance obtained from the reflectance image. For NIR, the simulated reflectance was slightly underestimated for all growth periods. The root mean square errors (RMSEs) between the obtained and the simulated reflectance values were approximately 10% to 20% of the averages.
Vegetation indices calculated using the simulated reflectance tended to be underestimated relative to those obtained from the reflectance image (Table 3). The RMSEs between the obtained and the simulated vegetation indices were the greatest for SR and the least for NDVI.

Discussion
In this study, FLiES, which is a tool based on a radiative transfer model and has been applied to crop canopies in previous studies [14,30], was used to simulate canopy reflectance for a paddy field for UAV monitoring. FLiES can simulate the process of reflection, transmission and absorption in a vegetation canopy and absorption by atmospheric particles in the air for each spectral band. A reasonable value of reflectance can be obtained by setting parameters properly in FLiES, including leaf and soil reflectance, LAD, plant height and coordinates of plants. In this study, the parameters were set based on the field measurements for paddy fields, and the effect of solar radiation conditions on the canopy reflectance and vegetation indices were evaluated by the simulation.
The canopy reflectance obtained from the reflectance image made from the aerial multispectral images taken by UAV for the measured paddy fields varied even at the points where LAI showed the same value. However, the simulated canopy reflectance will not vary if LAD and plant height do not change. The reflectance difference might be caused by a measurement error for LAI, which was estimated to be 30% in the paddy using the plant canopy analyzer LAI-2200 [36]. The fewer parameters set in the simulation may also explain this difference. The greater difference at lower LAI suggested that the canopy reflectance obtained from the reflectance image was influenced by the background conditions, such as dry soil, wet soil and water surface, because the influence of background on canopy reflectance would be greater than that of LAI at lower LAI. Modeling of spatial leaf distribution in the simulation might also not be accurate, even though the coordinates of plants, the distance between rows and the hill interval were considered. Evaluation of the background and spatial leaf distribution at each measured point is needed to precisely simulate the canopy reflectance. However, the average of the obtained reflectance and vegetation indices for each growth period were similar to the simulation results, and RMSE of NDVI decreased as NDVI became larger (Table 3), which agreed with the results in Borgogno-Mondino et al. [23]. These results suggest that the simulation is accurate enough to evaluate the average behavior of photons and the effect of solar radiation conditions on canopy reflectance.
The simulation with various solar radiation conditions conducted in this study suggests that two desirable conditions are needed for canopy reflectance measurement (see Figure 2). First, aerial multispectral images should be collected during cloudy weather when the proportion of diffuse light is high. Second, aerial multispectral images should be collected from 45 to 65 degrees of solar zenith angle during sunny weather when the proportion of diffuse light is low. Collection of aerial multispectral images should be avoided when the solar zenith angle is very large or small. These conditions are consistent with those suitable for the measurement of LAI using a plant canopy analyzer. In the case of a plant canopy analyzer, early morning is also suitable because the direction can be selected to avoid the effect of direct sunlight. However, because the collection of aerial multispectral images cannot select the measurement direction (i.e., nadir), early morning is not suitable for acquiring images. These limitations for the collection of aerial multispectral images decrease the flexibility benefit of UAV. Accordingly, the information needed to estimate the effects of solar radiation conditions, such as time, date, latitude and cloud coverage, should be recorded at the image collection. The simulation by FLiES suggests that the change in canopy reflectance caused by solar radiation conditions can be calibrated.
The calibration method will be tested in a future study.
The vegetation index is often utilized to evaluate crop growth. The simulation results in this study showed that NDVI is the most robust among the three vegetation indices tested across various solar radiation conditions. This may explain why NDVI is employed by many studies. However, NDVI saturates at the least LAI. SR and EVI2 hardly saturate when compared with NDVI, therefore, these indices would be better to evaluate crop growth if the change in canopy reflectance caused by solar radiation conditions was calibrated. The simulation results also indicate that SR is sensitive to solar radiation conditions. The sensitivity of SR is derived from Equation (1), of which the numerator and denominator directly enhance the variation of NIR and RED, respectively ( Figure 3). Consequently, EVI2 seems to be the most suitable index to evaluate crop growth among the three vegetation indices.

Conclusions
To clarify how solar radiation conditions effect canopy reflectance in aerial multispectral images obtained using UAV, a canopy reflectance simulation using a radiative transfer model was conducted for various solar radiation conditions and canopy structure parameters. The results showed that the canopy reflectance in both RED and NIR varied depending on the proportion of diffuse light for incident sunlight and the solar zenith angle. The variation in canopy reflectance was greater when the proportion of diffuse light was low (0.0 and 0.5) and the solar zenith angle was large (85 degrees) and small (5 and 25 degrees). The variation in canopy reflectance inevitably affected vegetation indices (SR, EVI2 and NDVI). SR was the most sensitive, and NDVI was the most robust for the solar radiation conditions. However, NDVI was saturated at the least LAI. These results suggest that EVI2 is the most suitable index to evaluate crop growth of the three vegetation indices. These simulation results were similar to the canopy reflectance and vegetation indices obtained from the aerial multispectral images collected by UAV for farmers' paddy fields.
Desirable conditions for the collection of aerial multispectral images, such as a high proportion of diffuse light for incident sunlight during cloudy weather, would decrease the variability caused by solar radiation conditions. In the case of undesirable conditions, a record of solar radiation conditions is required for the collection of aerial multispectral images. From the viewpoint of promoting the use of UAV in agricultural fields, the development of calibration methods to reduce the influence of solar radiation conditions is necessary to maintain the flexibility of UAV use. Acknowledgments: Thanks to Sendai Arahama, an agricultural producers' cooperative corporation, for providing the paddy fields for field measurements and aerial image collection. Thanks to all members of the Crop Science Laboratory of the Graduate School of Agricultural Science, Tohoku University.

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations
The following abbreviations are used in this manuscript. EVI2 Enhanced