A Robust Vegetation Index based on Different UAV RGB Images to Estimate SPAD Values of Naked Barley Leaves

: Chlorophyll content in plant leaves is an essential indicator of the growth condition and the fertilization management effect of naked barley crops. The soil plant analysis development (SPAD) values strongly correlate with leaf chlorophyll contents. Unmanned Aerial Vehicles (UAV) can provide an efficient way to retrieve SPAD values on a relatively large scale with a high temporal resolution. But the UAV mounted with high-cost multispectral or hyperspectral sensors may be a tremendous economic burden for smallholder farmers. To overcome this shortcoming, we investigated the potential of UAV mounted with a commercial digital camera for estimating the SPAD values of naked barley leaves. We related 21 color-based vegetation indices (VIs) calculated from UAV images acquired from two flight heights (6.0 m and 50.0 m above ground level) in four different growth stages with SPAD values. Our results indicated that vegetation extraction and naked barley ears mask could improve the correlation between image-calculated vegetation indices and SPAD values. The VIs of ‘L*,’ ‘b*,’ ‘G-B’ and ‘2G-R-B’ showed significant correlations with SPAD values of naked barley leaves at both flight heights. The validation of the regression model showed that the index of ‘G-B’ could be regarded as the most robust vegetation index for predicting the SPAD values of naked barley leaves for different images and different flight heights. Our study demonstrated that the UAV mounted with a commercial camera has great potentiality in retrieving SPAD values of naked barley leaves under unstable photography conditions. It is significant for farmers to take advantage of the cheap measurement system to monitor crops.


Introduction
Barley is the fourth largest cereal crop after maize, wheat and rice globally, with 7 percent of global cereal production. Barley production around the world is 159.0 million tons in 2019 [1]. It is a major food source for the people in the cool and semi-arid areas of the world, where wheat and other cereals are less adapted [2]. Due to the relatively harsh growing conditions, some suitable managements, especially nitrogen fertilizer, are necessary to sustain healthy growth and maximum yield production [3]. Moreover, leaf chlorophyll content has been proved to be related to nitrogen content in the crops and can be used as a fertilizer indicator during crop management [4][5][6][7][8]. The traditional chemical analysis of chlorophyll content has limitations of labor intensity, destruction and small measurement area. The soil and plant analysis development (SPAD) values measured by SPAD-502Plus Chlorophyll Meter (KONICA MINOLTA, Japan) can provide an accurate assessment of chlorophyll content in plant leaf in a non-destructive way [9][10][11][12][13]. But the measurement area of SPAD-502 is still limited and the meter can only be conducted on leaves one by one [14]. The method that could estimate SPAD values of leaves on a large scale and short time would help crop fertilizer management in the field.
Nowadays, due to the development in color image acquirement, the color-based vegetation indices (VIs) from digital color cameras have been widely used in diagnosing crop physiological characteristics, including texture and shape features, leaf nitrogen and chlorophyll content [14][15][16][17]. For example, the vegetation indices based on red, green and blue (RGB) color system, such as 'G/R,' 'R-B,' '2G-R-B' and '(R-B)/(R+B)' are usually called "about-greenness" indices for analyzing plant characteristics, for estimating leaf chlorophyll and nitrogen contents or for distinguishing plant from the non-plant background [18][19][20][21][22]. In addition to the RGB-based vegetation indices, the vegetation indices based on other color systems, like International Commission on Illumination (CIE) L*a*b* color system and hue-saturation-intensity color system, also have been used in estimating vegetation status, like diagnosing macronutrient contents (nitrogen, phosphorus, potassium and magnesium) in crop leaves [23,24]. For a digital color camera, the image appearance (such as brightness and color temperature) would be influenced easily by illumination. The stable light condition could generally provide optimal light intensity and stability for the digital image capture and then the digital image would be convenient for the data analysis. On the contrary, the sunlight intensity changes would lead to different image appearances in different images [25].
In recent years, Unmanned Aerial Vehicles (UAV) have received more attention for their capability of covering a relatively large area in a short time and unlimited capture time, which could minimize the measurement error caused by changes in environmental factors [26][27][28][29][30]. Due to the advantages of low cost, portable and high spatial and temporal resolution (compared to satellite remote sensing), UAV has been widely used in agriculture to estimate crop traits, such as chlorophyll and nitrogen content, plant height, above-ground biomass, yield and fractional vegetation cover [31][32][33][34][35][36]. However, UAV mounted with multispectral or hyperspectral sensors would lead to high expense, which is not friendly to the ordinary farmers. Hence, the low-cost UAV systems mounted with a commercial digital color camera would be an alternative and cost-prohibitive method to make UAV popular among average farmers [37][38][39].
A lot of previous studies about retrieving crop characteristics (like estimating plant growth parameters, biomass, nitrogen and chlorophyll content and so on) by UAV-derived images, which included red, green and blue channels, were mostly based on the image acquired at the same flight height [40][41][42][43][44][45]. The photography time of all UAV images in these studies was no more than 30 minutes to ensure the images could be successfully overlapped with others; thus, the photography condition (like sunlight intensity and solar elevation angle) of all images was similar. There was no need to consider these variances among different images during image analysis. But for the study based on UAV-images acquired on different days, the solar intensity and elevation angle among images may not be consistent and become an influencing factor on the biophysical parameter retrieval from the image because of the necessity of combining all images. Moreover, different flight heights would also affect analysis results based on UAV-images [46,47]. However, the different environments during the image acquirements, such as photography conditions of solar and flight heights, on the SPAD values estimation were not profoundly studied. Therefore, a color-based vegetation index that can estimate leaf SPAD values independent of the solar condition and the photography height during UAV image acquisition by a commercial digital color camera is necessary.
Therefore, the objectives of this study were (1) to establish the optimal image processing methods for deriving good relationships between naked barley leaf SPAD values and vegetation indices using different UAV images and (2) to find a robust vegetation index for estimating the naked barley leaf SPAD values using UAV images taken at different photography conditions of solar and fight heights.

Experimental Field Design
The study area is located at the experimental station of Ehime Research Institute of Agriculture, Forestry and Fisheries in Matsuyama City, Ehime Prefecture, Japan (33°58ʹ N, 132°48ʹ E, World Geodetic System 1984). The area is belonged to a temperate zone, with an average yearly temperature of 16 ℃ and the precipitation of 1300 mm. The experiment was conducted from December 2018 to March 2019.
The experiment included 6 experimental groups (Group A, Group B, Group C, Group D, Group E, Group F) for different fertilization treatments. The fertilization details of each group are shown in Table 1. Each group has an area of 5.4 m × 9.0 m and repeats three times (Repeat 1, Repeat 2, Repeat 3). The field map is shown in Figure 1.
The naked barley measured in the experiment was Haruhimeboshi (Hordeum vulgare L. var. nudum). The sowing rate was 10 kg/10a.

Data Collection and Analysis
The SPAD values of naked barley leaves and UAV images of the canopy were acquired at different growth stages. The timings of collecting SPAD values and capturing UAV images and other managements (fertilizer application, sowing and harvesting) are described in Figure 2. The data was acquired at the growth stage (December), stem elongation stage (January and February) and flowering stage (March) according to the BBCH (Biologische Bundesanstalt, Bundessortenamt und CHemische Industrie)-scale [48]. SPAD-502Plus Chlorophyll Meter was used to measure the SPAD values of naked barley leaves in the field. The meter has a measuring area of 2 mm × 3 mm and a measuring accuracy within 1.0 SPAD unit. The average SPAD value by 3 random measurements at different places in one leaf was used as the leaf's measured SPAD value. The average SPAD value by 10 naked barley leaves in each group was used as the experimental groupʹs measured SPAD value.
DJI Phantom3 Professional (Da-Jiang Innovations Science and Technology Co., China) was the UAV system to acquire naked barley canopy images in the RGB wavelengths in this study. The camera uploaded on the UAV uses a 1/2.3 complementary metal-oxide-semiconductor (CMOS) with 12.4 megapixels (4000 × 3000) and the shooting direction is perpendicular to it. The field of view of the camera is 94˚, the focal length is 20mm, the f-number is f/2.8 and focus at the infinite point (∞). The camera uses automatic exposure during the image acquirement. The reflectance was not obtained by the reference normalization because the atmospheric effect was regarded very little due to the low flight heights in this study. Thus, the DN values of R, G and B channels of UAV images represented the relative radiance. The images were obtained at two different flight heights above ground level of 6.0 m and 50.0 m with an accuracy of flight altitude of 0.1 m. For the images at 6.0 m, each image only included one experimental group and there were no overlapped parts among different group images. Whereas the images at 50.0 m included all experimental groups (18 groups) for one day.
The measured SPAD value in the field and UAV images at 6.0 m and 50.0 m of one group were used as one set of data samples in the study. There were 18 samples per day and 4 flight days. Nevertheless, the data of 'Repeat 2' of 'Group C' on 8th December 2018 was deleted from the sample set because of the damaged UAV image at 6.0 m, so that the total sample number was 71.

Flowchart of Data Processing
The whole data processing was illustrated in Figure 3. Firstly, the image classification by Maximum likelihood classification (MLC) was conducted to extract vegetation parts on the UAV images of both flight heights; secondly, the naked barley ear was subsequently masked for the images at 6.0 m; thirdly, the VIs were calculated for the original, vegetation extracted without the ear masked and vegetation extracted with the ear masked images; fourthly, linear regression models were developed to relate the VIs with the field measured SPAD values; fifthly, the SPAD values prediction models were validated by an independent dataset; lastly, a robust VI for predicting SPAD values was selected for different photography conditions (including solar and flight height).

Image Segmentation
MLC was used to extract the vegetation cover by classifying UAV-derived RGB images at both flight heights by the ENVI 5.5 (Harris Geospatial Solutions) software [49]. At first, for both flight heights, the UAV-images were classified into three categories: naked barley, soil and shadow. Then, the soil and shadow parts were masked and the vegetation cover was left in the images. Moreover, the naked barley ear was additionally masked after vegetation extraction for the images of 6.

Color-Based Vegetation Indices
A 5 m × 4.5 m image area was used for VIs calculation for each group. Due to the different photography area in one image at different flight heights, 3.6 megapixels (2000 ×1800) were used for each experimental group at 6.0 m flight height and 28.8 thousand pixels (180 × 160) were used for each group at 50.0 m flight height. The marginal areas of each group were cut from the images to eliminate the effect of fertilization on naked barley growth in different groups because the fertilization details were different among the neighboring groups. The cut-images would be called UAV images below.
The CIE L*a*b* color system is designed to approximate human vision. It describes color by three values: (1) L* represents the human perception of lightness from black (0) to white (100), (2) a* represents the visual perception from green (-a*) to red (+a*) and (3) b* means for the visual perception from blue (-b*) to yellow (+b*). L*, a* and b* are computed from tristimulus values X, Y and Z as following equations: (1) where , and describe a specified gray object-color stimulus [23,24,57].  48 samples of 'Repeat 1' and 'Repeat 3' were used to detect the relationship between VIs and measured SPAD values by Pearson correlation and regression analyses after checking the normality assumption. Linear regression models were fitted to each VI. The significance of linear regressions was evaluated using Student's t-test at 95% confidence levels. The data analysis was done by IBM SPSS Statistics (International Business Machines Corporation). 23 samples of 'Repeat 2' were used to validate regression models.
For each regression model, the root mean square error (RMSE) and normalized root mean square error (NRMSE) were used to validate using 23 samples of 'Repeat 2.' The equations of RMSE and NRMSE are presented in Equations (8) -(9).
in which, is the measured SPAD value of sample , ŷ is the predicted SPAD value of sample and is the number of validation samples.

SPAD Values for Different Fertilizer Treatments
The measured SPAD values of each experimental group in the field at different dates were shown in Figure 5. The lowest SPAD value (26.30) occurred at Repeat 2, Group C, January 2019 and the highest (55.95) appeared at Repeat 1, Group C, March 2019. The SPAD value range was not much different from that in other researches, for example, 2~60 for sugar beet and 2.88~53.43 for potato [58,59]. For the SPAD values measured on the same date, the maximum difference occurred on 8th January 2019, from 26.30 (Repeat 2, Group C) to 36.30 (Repeat 1, Group B). The results of one-way analysis of variance of SPAD values for different fertilizer treatments and growth stages were shown in Table 3. Only the SPAD values in Jan. 2019 showed a significant difference for different fertilizer treatments. In contrast, there were significant differences among different growth stages for all fertilizer treatments. The overview for all experimental groups is displayed in Figure 6 by the UAV-RGB images derived at the flight height of 50.0 m. The parts corresponding to the low SPAD values showed light green, even close to yellow, while the regions with high SPAD values displayed green or dark green in the images ( Figure 6). The maximum difference in the scene due to the fertilizer treatments occurred on 8th January 2019 ( Figure 6(b)).
However, there are some slight color differences among the same experimental groups on the same date or even in the same sample plot. Generally, the solid fertilizer used in our study can only be absorbed by the crop after being dissolved. Precipitation is one of the most critical impact factors on the degree of fertilizer dissolved under the rainwater irrigation system. Moreover, the naked barley growth density for different experimental groups in our study was not the same, caused by the inhomogeneity of nutrient, moisture and other soil characteristics. All these factors would lead to different fertilizer utilization rates of the naked barley, even for the same fertilizer treatment, in the field. Therefore, growth status variances occurred to the same fertilizer treatment experimental groups in our study and showed different greenness in the overview image on the same date.

Optimal Image Processing for Best Relating SPAD Values with Vegetation Indices
The VIs were calculated from the original image, the image with vegetation extraction without naked barley ears mask and the image with the vegetation extraction and naked barley ears mask, respectively, for the 6.0 m height images. Simultaneously, the vegetation indices were calculated only on the original image and the image with vegetation extraction for the 50.0 m height because the ears mask was not applied on it. Correlation between measured SPAD values and VIs of sample sets of 'Repeat 1' and 'Repeat 3' were displayed by correlation coefficient (r) in Table 4. The correlation coefficients after applying vegetation extraction were higher than that in the original imagery data for most VIs at both flight heights. Moreover, for the VIs calculated from UAV images at 6.0 m, the correlation coefficients showed even higher on the images after vegetation extraction and naked barley ears masked than that with the only vegetation extraction. In other words, the VIs had a better correlation with SPAD values of naked barley leaves after applying both vegetation extraction and naked barley ears masked.

Relationships between SPAD Values and Vegetation Indices
The different VIs showed various relationships with SPAD values. It seemed that most of the VIs did not perform consistently on the images from different flight heights.  Figure 7 and 8, respectively. Among these 4 VIs in Figure 7 and Figure 8, the VI of 'G-B' performed the best correlation with SPAD values of naked barley leaves at both flight heights.

Validation of the Linear Regression Models for Predicting SPAD Values
The validation of the linear regression models for predicting SPAD values was conducted by an independent dataset of 'Repeat 2' (n=23). The RMSE and NRMSE of regression models were shown in Table 5. The table suggested that the RMSE and NRMSE of the indices, which performed better correlations with SPAD values (Table 4) Figure 11 with Figure 6, we could find that the SPAD value image roughly corresponded to the colors and tones of the RGB image. The pixels with low SPAD values in Figure 11 corresponded to the light green area, even similar to yellow color, in Figure 6. The pixels with high SPAD values in Figure 11 corresponded to the green or dark green area in Figure 6.

The SPAD Value Estimation Improvements by Vegetation Extraction and Naked Barley Ears Mask
The vegetation extraction and naked barley ears mask have improved the SPAD value estimation (Table 4) because they removed most of the soil, shadow and naked barley ears in the image. The soil, shadow and naked barley ears parts were considered as the background noise in our study because that they were irrelevant to SPAD values of naked barley and would affect the vegetation discrimination in the image [60,61]. Thus the image-based VIs can have a stronger correlation with SPAD values after the background was masked. Therefore, we consider that the vegetation extraction (soil and shadow masked) and naked barley ears mask as an optimal image processing method for correlating color-based vegetation indices calculated from UAV images with SPAD values of naked barley leaves.

Performance of the Color-Based VIs on Estimating SPAD Values
According to Table 4, we could find that very high differences occurred on the rela-  (Figure 6), which would lead to some errors in the image-based vegetation indices calculation on the images at 50.0 m height [62][63][64]. This inconsistent brightness in the whole image may be caused by the unstable flight status and camera's internal parameter adjustment in a large photography area of the commercial digital camera mounted on the UAV in our study. The UAV images at 50.0 m in our study had a large image area, including the experimental area with low reflectance and some bare sandy soil with a high reflectance located in the north of our experimental area. In the south area of the image at 50.0 m, only farmland with low reflectance was included. For the north area in the image, the camera's internal parameters might be difficult to adjust due to the clear difference in the brightness between the bare sandy soil and farmland in the north of the image, especially on clear days. Therefore, the inconsistent brightness occurred to the image at 50.0 m due to the digital camera's internal parameter adjustment during photographing. However, for the images at 6.0 m, there was only the experimental area in each image. For the UAV images captured at the flight height of 6.0 m, compared to the images at 50.0 m, a smaller photographic area was performed in each image. There was only the experimental area in each image of 6.0 m. It may be easy for the commercial digital camera to adjust its internal parameters to a more suitable condition to acquire the image at a lower flight height. As a result, there was no apparent non-uniform brightness in the images at 6.0 m. Therefore, the distinctly inconsistent brightness throughout the whole RGB image at the flight height of 50.0 m could be considered as the reason for variances between the relationship of SPAD values and the index of 'R,' 'G,' 'B' in two different flight heights in our study.
For the color-based VIs calculated from UAV images at 6.0 m and 50.0 m, the index of 'L*,' 'b*,' 'G-B' and '2G-R-B' performed significant correlations with SPAD values of naked barley leaves at these two flight heights in our study and 'G-B' was the most robust one among them. For the VIs of 'G-B' and '2G-R-B,' they could be regarded as corrected values for the form of subtracting the 'B' channel value from the 'G' channel value or subtracting the 'R' and 'B' channel values from the 'G' channel value during the VI calculations. They may eliminate the brightness unbalance within one image for different experimental groups at 50.0 m height in our study. Consequently, the correlations between the index of 'G-B' and '2G-R-B' with SPAD values were strong at 50.0 m, even under the whole imageʹs different brightness. Further research that focuses on the brightness calibration among whole UAV-derived images should be expected.
The index of 'L*' and 'b*,' which were based on CIE L*a*b* color system, also showed reasonable performance for estimating SPAD values. During the calculation from RGB values (Equations 1-7), the 'G' value took up a large proportion of the index of 'L*' and the index of 'b*' responded to the index of 'G-B.' Also, the index of 'L*' and 'b*' was divided by a gray board during calculating, which might normalize the index. Therefore, the index of 'L*' and 'b*' also performed significant correlations with SPAD values.
Moreover, for the relationship of 'R,' 'G,' 'B' indices with SPAD values of naked barley leaves in our study, there was also one more point that should be noted. The order of relevancy for these three single channels with SPAD values was not the same in our studyʹs two flight heights. The 'R' and 'G' had a much higher correlation with SPAD values than 'B' and 'R' was a little higher than 'G' for the UAV-image at the flight height of 6.0 m, whereas 'G' had much higher correlation than 'R' and 'B' and 'R' was a few higher than 'B' for the UAV-image at 50.0 m. Furthermore, neither of these two correlations was consistent with a previous study about correlation analysis of total chlorophyll content with 'R,' 'G' and 'B' indices, which were derived from the multispectral TV camera and conducted under constant lighting in the laboratory through various interference filters in a silicon vidicon camera [65]. The index of 'B' showed the best correlation and 'R' was the third in the previous study. Hunt has indicated that the vegetation indices, which usually showed good performance at the leaf scale, may not perform well at the canopy level [54]. Generally, the 'B' value of targets, one of the RGB color system parameters, seldom can be stably measured under solar light. The 'B' value is a channel easily affected by the atmospheric conditions because of the scattering; meanwhile, the 'R' and 'G' are not easily affected by atmospheric scattering when compared with the 'B' channel, especially on a clear day [66]. The indices used in our study were based on the UAV-derived images in the field under the sunlight. The reflected light from the naked barley canopy and the scattered light from the atmosphere were included during vegetation index calculation for the 'B' channelʹs value. This would lower the precision of 'B' values in our images and lead to the instability of chlorophyll content estimation accuracy among different light sources during photography.
Although the commercial digital color camera has the advantages of low cost, convenient operation and simple data processing, exposure time and internal parameters adjustment, which would be set based on the weather conditions, would significantly affect image quality such as brightness and color temperature [67]. The unstable photographic condition, including weather and atmospheric conditions described above, may result in the instability of a single channel's color-based index calculation. On the contrary, the vegetation indices of 'G-B' and '2G-R-B,' the adjusted indices, obtained a relatively high correlation with SPAD values of naked barley leaves in our study.

Conclusions
The present study suggested that the application of vegetation extraction and naked barley ears mask can improve the correlation between color-based vegetation indices calculated from UAV-mounted commercial digital camera and SPAD values of the naked barley leaves. Although vegetation indices had different correlations between SPAD values for different UAV flight heights and photography conditions, 'L*,' 'b*,' 'G-B' and '2G-R-B' had significant correlations with SPAD values under different conditions. The index of 'G-B' can be regarded as the best vegetation index for estimating SPAD values of naked barley leaf from UAV-derived RGB images in our study because of the excellent performance in the validation. This vegetation index could be used for different photography conditions (including sunlight condition, camera's internal parameter and UAV flight height).
This research provided a method to estimate SPAD values of naked barley in a wide range by UAV. The method could reduce the influence of UAV photography conditions, which can make environmental conditions for acquiring UAV images less strict and did not require multispectral or hyperspectral sensors, which are very high expenses for ordinary farmers. The commercial digital camera usage, which showed a reasonable estimation of SPAD values of the naked barley crop, could be regarded as a low-cost method. Hence, its affordability may increase its application in the crop monitoring practice.
From a longer-term perspective, more deeply research about the relationship between crop biophysical parameters (not only SPAD values but also above-ground biomass, yield and so on) and color-based vegetation indices derived from UAV-mounted commercial digital camera and crop biophysical parameters is expected. Funding: This research received no external funding.

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