Calibration and Image Processing of Aerial Thermal Image for UAV Application in Crop Water Stress Estimation

ISAE-Supaéro, Université de Toulouse, 10 Avenue Edouard Belin, 31055 Toulouse, France Department of Biosystems Engineering, Seoul National University, 1 Gwanak-ro, Gwanak-gu, Seoul 08826, Republic of Korea Global Smart Farm Convergence Major, Seoul National University, 1 Gwanak-ro, Gwanak-gu, Seoul 08826, Republic of Korea Research Institute of Agriculture and Life Sciences, Seoul National University, 1 Gwanak-ro, Gwanak-gu, Seoul 08826, Republic of Korea


Introduction
The era of smart farms and precision agriculture has arrived, where various technologies are being applied, including artificial intelligence, Internet of Things, information and communication technologies, autonomous vehicles, and unmanned aerial vehicles (UAVs) [1,2]. Among these rapidly developing technologies, UAVs have numerous advantages and are thus used in various fields. UAVs have multiple advantages over other remote sensing airborne platforms and satellites, including data capture with a higher resolution, lower cost, and higher flexibility [1, [3][4][5]. In addition, the use of UAVs provides numerous advantages over ground-based methods, which tend to be more time consuming, costly, and often destructive [6].
Among the above UAV technologies, drones with various remote sensors are being widely used for remote sensing of crops [7]. For example, UAV remote sensing technologies are being used for crop monitoring, irrigation and fertilization management, pathogen detection, yield estimation, growth vigor assessment, and water stress detection with the use of drones and infrared thermal cameras [8]. The use of UAVs with cameras and sensors can be an alternative to both ground-based measurements and airborne/satellite data measurements, which are more time consuming [9]. The use of UAV remote sensing enables a nondestructive and rapid data capture, compared to ground-based measurements. Furthermore, the UAV remote sensing technology provides a competitive advantage in terms of cost and turnaround time over airborne and satellite solutions [10]. In addition, the drone can fly safely at low altitudes and capture highspatial-resolution images at a lower price, compared to even manned aircrafts such as helicopters [11]. Another advantage of using a UAV is the flexibility of its applications, which enables multitemporal studies by a small revisiting time [12] as well as higher independence of climatic variables. Owing to the multiple advantages of UAVs over other remote sensing platforms, the UAV remote sensing technology could be widely used in precision agriculture, where the collection, processing, and analysis of actionable data are very important.
The evaluation of water stress and irrigation management is considered one of the most important applications in precision agriculture [13]. The ability to assess the water stress status of plants is crucial, as it allows the implementation of site-specific irrigation systems that could supply the optimal amount of water for the actual need of the plants [14]. Moreover, a higher degree of agricultural automation can be achieved [15]. The plant water stress evaluation methods include indicators based on measurements of both soil and plants [16]. Regarding plant-based measures for water stress, nondestructive methods have been proposed, which include the infrared thermal imaging technology with both ground-based methods and UAV-based methods. For example, UAVs equipped with thermal cameras have been utilized to assess the agricultural crop's response to water stress [17]. Consequently, UAV-based aerial thermal imaging can be a very valuable method for the assessment of the water status of plants [18]. Regarding the evaluation of water stress response in plants and trees, traditional ground-based techniques often require measurement of stomatal conductance, involving direct contact with the leaves, which is time consuming [6]. However, with the use of UAVs and thermal cameras, it is possible to capture large-scale aerial thermal images of plants and then utilize the relationship between the canopy temperature and stomatal conductance [10].
To assess the validity and feasibility of aerial thermal imaging, as an approach to assess the water status in plants, numerous studies involved water treatment experiments to analyze the crop response with the use of UAVs and thermal cameras. For example, Gomez-Candon et al. [19] assessed the water stress response in an apple tree field using remotely sensed thermal images captured with a drone, where half of the trees were well irrigated, while the other half were subjected to progressive summer water stress. Ludovisi et al. [18] utilized the UAV-based remote thermal imaging technology for field phenotyping of black poplar response to water stress. In this case, aerial thermal imaging has been employed because of the efficient and high-resolution (leaf level) nondestructive monitoring of the genotype performance in large populations. Grant et al. [17] analyzed whether significant differences exist in water tension, leaf stomatal conductance, and thermal index according to differ-ences in irrigation conditions for grapevines, soybeans, and flowers and confirmed the linearity between the indices. De Jonge et al. [20] quantified the crop water stress using six types of water stress index (crop water stress index (CWSI), degrees above nonstressed (DANS), degrees above canopy threshold (DACT), time temperature threshold (TTT), integrated DANS (IDANS), and integrated DACT (IDACT)). They confirmed that the early afternoon time was the best time for the calculation of the water stress index and that all water stress indices have linear relationships with the soil moisture at high temperatures.
Although multiple advantages exist over aerial infrared thermography, the raw aerial thermal images need to be calibrated because of the difference between the measured temperature by the aerial thermal images and actual temperature of the target, owing to several factors, including the altitude of the UAV and atmospheric conditions such as the humidity and solar irradiation [21,22]. Furthermore, the aerial thermal images need to be postprocessed because they contain noise components such as data related to the soil temperature measurement and other objects.
In this study, we propose an algorithm for the calibration of aerial thermal images regardless of altitude and atmospheric conditions and thermal image processing methods, which enable to selectively extract only the temperature of the canopy from the aerial thermal images, while excluding other irrelevant elements such as the soil and objects. We then calculated the DANS water stress index using the canopy temperature extracted from the aerial thermal images subjected to calibration and processing and analyzed the relationships between the irrigation conditions and calculated DANS water stress index.

Material and Methods
2.1. Calibration of the Aerial Thermal Image. For the calibration of the aerial thermal image, we used a blackbody with a uniform surface temperature. The blackbody was composed of (1) a radiator fabricated using an aluminum plate, (2) carbon-fiber-based meshed heat pad, (3) insulating material fabricated using a polyethylene foam, and (4) outer frame of an acrylic material. The width, length, and height of the blackbody system were 730, 730, and 50 mm, respectively, while its mass was 5 kg, which made it portable [23]. For the measurement of aerial thermal images of the blackbody, a drone (DJI Phantom 3, DJI, China) equipped with a thermal camera (FLIR VUE PRO R, FLIR Systems Inc., USA) was utilized, as shown in Figure 1. Even though this thermal camera has some limitations of calibration shift and vignetting during its operation, it was corrected by recalibration with a shutter and side overlapping of 90% during the aerial shooting.
The calibration of the aerial thermal image by a blackbody was performed as follows. (1) At a constant radiant heat from the blackbody on the ground, the UAV equipped with the thermal camera captured thermal images of the blackbody every 10 m, from 10 up to 80 m, as shown in Figure 2.
(2) The blackbody's temperature was measured on the ground, using resistance temperature detector (RTD) 2 Journal of Sensors contact-type temperature sensors in contact with the blackbody radiator during the entire experiment. (3) During the data capture of the aerial thermal images, the atmospheric temperature, relative humidity, and solar irradiance were measured on the ground. (4) The aerial thermal images and other measured variables were then used for calibration using the Gaussian process regression (GPR) method. The regression was based on Eq. (1), which predicts the difference between the blackbody temperature measured by the RTD temperature sensors (TbG) and blackbody temper-ature measured by the UAV (TbA) with independent variables of altitude (Alt), atmospheric temperature (Temp), relative humidity (Humd), and solar irradiance (Irrad), TbG-TbA = fun Alt, Temp, Humd, Irrad ð Þ : To proceed with the regression analysis, the independence between the independent variables needs to be confirmed with the preceding conditions. In this study, we confirmed the independence between independent variables by analyzing the scatter matrix plot and variation inflation factor (VIF). When the linearity between the independent variables in the scatter matrix plot is observed or strong linearity is estimated during the VIF calculation, we need to remove the independent variable that affects the dependency between the independent variables. If the VIF is larger than 10, the multicollinearity is high. If the VIF is higher than 5, it is not that high, but the use of this variable would require a careful analysis [24].
The GPR method was applied to the collected data to create a temperature correction model for the altitude and environmental variables (atmospheric temperature, relative humidity, and solar irradiance). The GPR method provides a nonparametric approach to predict a model from data without assuming the structure of the model, rather than a parametric approach that predetermines the model and predicts a variable as a linear regression equation, which may be obtained as a statistical expression, as shown by Eqs. (2)-(5). It is a method of predicting the mean and covariance of the normal distribution where the dependent variable y * is located when a new independent variable X * is given based on the independent variable X and dependent variable y of the existing training  3 Journal of Sensors data [25]. The conditional distribution of y * can be expressed as a multivariate normal distribution, where μ * is the mean, Σ * is the covariance matrix, σ 2 is the assumed noise level, and I is the identity matrix [26].
The rational quadratic kernel, K, is applied to the calculation of the covariance matrix, where three hyperparameters, σ f 2 (overall variance), α (mixture scale), and l (length scale), were adjusted with the given dataset [27].
The total number of training and test sets was 2,912, while the percentages of the observations in the training   Journal of Sensors and test sets were set to 80% and 20%, respectively. For each set, we evaluated R 2 (coefficient of determination) and root mean square error (RMSE) and verified the regression model by analyzing the following characteristics of the residuals: linearity, homoscedasticity, and independence. The linearity and homoscedasticity of the residuals were evaluated using    [28]. The GPR model and its evaluation methods were developed using a commercial programming tool (MATLAB ver. R2016, MathWorks Inc., Natick, MA, USA). Figure 3, aerial thermal images of crops were collected at the experimental farm located in National Institute of Horticultural and Herbal Science (Rural Development Administration, Korea). Aerial thermal image processing methods were realized to precisely extract the canopy temperature of fruit trees from aerial thermal images obtained by the UAV. The subjects of the experiment were 21 apple trees, located within an apple orchard (length: 40 m, width: 3 m). Number plates were placed at both sides of each tree to distinguish a target tree from all trees. The 21 apple trees were given water in different doses based on different percentages of potential evapotranspiration (ET) as follows: four trees under 100% ET, five trees under 75% ET, five trees under 50% ET, four trees not subjected to irrigation, and three trees not subjected to irrigation and rainfall. Aerial thermal images were collected four times in the morning and two times in the afternoon, which provided approximately 200 thermal images. The UAV was flown at an altitude of 10 m, with a side overlap of 90%. To compare the accuracies of the aerial thermal image processing methods, the leaf temperatures of every tree were measured manually by imposing a size of 20 × 20 pixel box on every tree in the collected aerial thermal images, as shown in Figure 3.

Aerial Thermal Image Processing. As shown in
We applied three types of aerial thermal image processing methods (Gaussian mixture model, Otsu algorithm, and Otsu algorithm after Gaussian blurring) to extract and classify the canopy temperature. As shown in Figure 4, in the Gaussian mixture model, n Gaussian distributions are mixed, which was employed to separate the unwanted background including soils from the measured aerial thermal images [29]. The Gaussian mixture model algorithm was applied under the following assumptions. First, the temperature of the canopy and background (soils) obeys a Gaussian distribution. Second, a mixture of two Gaussian distributions produces two different temperature distributions for the aerial thermal image. Third, the lower maximum of the two distri-butions was assumed to be the average temperature of the canopy, because the canopy temperature is generally lower than the soil temperature during the season of this experiment. Under these three assumptions, the average temperature of the canopy was extracted by the temperature distribution of the thermal image of each partition as the sum of two Gaussian distributions.
Binarization is a method of separating a subject into two parts using a threshold value. The temperature information separated by a section was regenerated into gray-scale images and then separated into canopy and background using the binarization method. The threshold value was determined using the Otsu algorithm. The Otsu algorithm is the most representative binarization method for determination of the threshold value. It can be used for the segmentation of thermal images [30]. The Otsu algorithm uses statistical values to determine the thresholds. It is an algorithm used to find gray values in the range of 0 to 255, which minimizes the cost function, where p1 and p2 represent the appearance probabilities of pixels corresponding to the first and second classes in the image, while V1 and V2 are the variances of gray values for the pixels of the first and second classes, respectively. As the Otsu algorithm is used to maximize the difference between the average values of the two classes, the iteration of the algorithm increases the gray value from 0 to 255 and changes p1 and p2 to find the gray value when the cost function is minimal [31]. In addition, we compared the Otsu algorithm and Otsu algorithm after Gaussian blurring. Gaussian blurring is a type of filtering method that removes noise generated by a Gaussian distribution. As it is assumed that the canopy temperature and noise obey Gaussian distributions, a significant difference in the temperature distribution exists when the noise generated by the Gaussian distribution is removed. Average temperature of the target crops as well as the standard deviations of temperature was analyzed and compared for the above three methods. For the aerial thermal image processing methods used in this study, we used a commercial software (FLIR tools, FLIR Ltd., USA), and we also developed a software using a commercial programming tool (MATLAB ver. R2016, MathWorks Inc., Natick, MA, USA).

DANS Index.
As mentioned in Introduction, several types of water stress index have been used, such as CWSI, DANS, DACT, TTT, IDANS, and IDACT [20]. Among them, the CWSI is the most popular water stress index because it is a normalized index including the environmental effects during its calculation. However, to calculate the CWSI, we need environmental information such as the atmospheric relative humidity in addition to the canopy  Journal of Sensors temperature. The information on environmental conditions cannot be easily acquired for most orchard farms or field cultivation farms [32,33]. Therefore, the DANS canopy index was employed in this study to evaluate the water stress based only on the canopy temperature. For most real farms, the measurement of only the canopy temperature, in addition to the relative humidity, can be very convenient. The DANS index shows the difference in the water stress at the same time without considering the environmental effect. It is calculated as the difference between the canopy temperature of the fruit tree exposed to water stress and canopy temperature of fruit trees not exposed to water stress, where T c is the leaf temperature at given time h and T cNS is the leaf temperature in crops that do not receive water stress at the same time [34,35]. After the DANS water stress index was calculated using the temperature correction method and aerial thermal image processing method, it is presented for different irrigation conditions to show the feasibility of the temperature correction and aerial thermal image processing methods proposed in this study.

Results and Discussion
3.1. Selection of the Variables. Before the GPR analysis, it is necessary to select the valid variables for the regression analysis among the input variables. This was performed by  Journal of Sensors checking the multicollinearity between the variables in two manners, through the scatter matrix and VIF. The multicollinearity was checked for all four measured variables (altitude, solar irradiance, atmospheric temperature, and relative humidity). In Figure 5, the relationships between variables are shown through the scatter plot matrix between variables, and Figure 6 shows a high linearity between the atmospheric temperature and relative humidity.
In addition, the multicollinearity of each variable was numerically determined by the VIF, a measure of the correlation between independent variables. The results are shown in Table 1. As shown in Figure 6 and Table 1 (VIF: 8.411), the atmospheric temperature can be regarded as dependent on other variables owing to the value larger than 5. Therefore, we performed a regression analysis excluding the atmospheric temperature variable. To confirm the independence between the remaining three variables, the VIFs of the variables except for the atmospheric temperature were calculated again. The results are shown in Table 2. Therefore, multicollinearity between the three variables (altitude, solar irradiance, relative humidity) was not observed according to the VIFs. The GPR analysis in this study was performed with these three variables, judged to be independent.

GPR Analysis
Results. Using 80% of the acquired data as a training set, the GPR training was carried out through a quadratic rational kernel. The regression equation was evaluated using the remaining 20% test set. In addition, the regression equation was evaluated by analyzing the residuals. The R 2 values of the regression analysis for the training and test sets were 0.9798 and 0.9685, respectively. The RMSEs for the training and test sets were 0.6680 and 0.8417°C, respectively. Figure 7 shows the predicted responses against the actual responses for the training and test sets. The RMSE of the test set was 0.8417°C, significantly lower than the error of 4.78°C, calculated before the GPR calibration. Thus, the GPR calibration was proven to be fairly effective. In addition, this method had a high R 2 for the training and test sets.
The residual analysis for this regression method was carried out for the following three properties of the residuals: linearity, homoscedasticity, and independence. The linearity and homoscedasticity of the residuals could be confirmed through Figure 8, showing graphs of the residuals for each predicted value of the training and test sets. As shown in Figure 8, except for the outliers in some sections, the training set is distributed evenly around zero of the y-axis and appears as a random noise. Therefore, the linearity of the residual could be confirmed. In addition, the homoscedasticity could be confirmed by the relatively uniform variance in all intervals. In the case of the test set, the outliers seem to have a large effect on the graph because of the lack of data. However, we confirmed that the test set was not significantly different from the training set.
To evaluate the independence between the residuals, a DW test [30] was performed, which yielded a value of 1.9926 for the training set and 1.8199 for the test set. As the value was close to 2, autocorrelation between the residuals was not observed. However, in this section, the sharing of the data for the training and test sets may have had an impact. Therefore, we proceeded with the learning for the entire dataset without separating the training and test sets. The DW test in this case yielded a value of 1.2814. Thus, autocorrelation existed, compared to the case of separating the training and test sets. However, compared to other regression methods, autocorrelation was not observed. Finally, the aerial thermal images were calibrated using the proposed GPR model. Figure 9 shows a representative thermal image corrected by the GPR model. The variables used for the GPR calibration of this thermal image were the shooting altitude, solar irradiance, relative humidity, and atmospheric temperature, and the values of them were 10 m, 958.33 W/m 2 , 45.7%, and 35.9°C, respectively.

Canopy Temperature for Image Processing Methods.
Three types of aerial thermal image processing method were applied to precisely extract the canopy temperature of the   Journal of Sensors fruit trees. These image processes were carried out after we cropped the tree of interest from the entire aerial thermal images, as shown in Figure 10. Figure 10(a) shows an aerial thermal image without the application of an image processing technique. Figures 10(b)-10(d) show the aerial thermal images processed by the Gaussian mixture model, Otsu algorithm, and Otsu binarization with Gaussian blurring, respectively. In these aerial thermal images, only the canopy temperature is presented, while other areas are reconstructed with zero value, shown as black.
The 21 apple trees were differently irrigated: four trees (Nos. 1-4) under 100% ET, five trees (Nos. 5-9) under 75% ET, five trees (Nos. [10][11][12][13][14] under 50% ET, four trees (Nos. [15][16][17][18] not subjected to irrigation, and last three trees (Nos. [19][20][21] not subjected to irrigation and rainfall. To evaluate the accuracies of the three aerial thermal image processing methods applied to extract the canopy temperature, the average and standard deviation of the canopy temperature extracted by each method were analyzed. The average value of the manually extracted canopy temperature for the aerial  Journal of Sensors thermal image without applying the image processing technique was 45.83°C, while the average values of the extracted canopy temperatures after applying the Gaussian mixture model, Otsu binarization method, and Otsu binarization with Gaussian blurring were 46.56, 46.97, and 46.88°C, respectively. The standard deviation averages were 1.25, 1.56, and 1.46, respectively. When the Gaussian mixture model method was applied, the standard deviation of the extracted canopy temperature was smallest and most similar to the manually analyzed canopy temperature. However, no significant differences existed between the image processing methods. Figure 11 shows the representative aerial thermal images of target crop for different five irrigation conditions ((a) under 100% ET, (b) under 75% ET, (c) under 50% ET, (d) not subjected to irrigation, (e) not subjected to irrigation and rainfall) on August 7th, 13th, 22th, and 29th. As shown in Figure 11, it was qualitatively observed that the canopy temperature of the crop under insufficient water supply (Figures 11(d) and 11(e)) is higher than that of the crop under sufficient irrigation (Figures 11(a) and 11(b)). Therefore, we employed three different thermal image processing methods (Gaussian mixture model, Otsu binarization method, Otsu binarization with Gaussian blurring), and we

11
Journal of Sensors quantitatively extracted the canopy temperature of the crops to make it easy and clear for the further data processing.  27,29,27, and 23°C, respectively. According to these temperature analyses, similar tendencies were observed for the measured canopy temperature and atmospheric daytime temperature, even though the differences between the canopy temperature and atmospheric daytime temperature were fairly larger than expected.

DANS Index for Different Irrigation Conditions.
In Figure 12, it is shown that the DANS index calculated by the canopy temperatures of target crops exposed to different irrigation conditions on different days. As shown in Figure 12, the DANS index increased according to the irrigation conditions (G1: under 100% ET; G2: under 75% ET; G3: under 50% ET; G4: no irrigation; G5: neither irrigation nor  Journal of Sensors rainfall). The DANS index was increased when the crops were exposed to a low irrigation on August 7th, 13th, and 22th. In addition, on August 29th, when the atmospheric temperature was lower than those on the other days, the water stress index did not exhibit a clear trend with respect to the irrigation conditions compared to the other days of experiments. Therefore, it was observed that the DANS index can indicate water stress conditions of crops using aerial thermal images, except when the atmospheric temperature abruptly decreases during the experimental period or on a cloudy day.

Conclusions
In this study, we realized an aerial thermal image calibration method and processing techniques to calibrate and extract the canopy temperature, which was then used to analyze the water stress level of fruit trees under different irrigation conditions. The calibration model was constructed to calibrate the aerial thermal images according to altitude and atmospheric conditions to match the actual ground temperature. The image processing technique was necessary to selectively extract only the temperature of the canopy from the aerial thermal images, while excluding irrelevant elements such as the soil and other objects. The image processing techniques included the Gaussian mixture model, Otsu binarization method, and Otsu binarization method after Gaussian blurring. The Gaussian mixture model provided the highest accuracy and stable results in extracting the canopy temperature. After the aerial thermal images were subjected to calibration and image processing, the DANS index was calculated using the canopy temperature of fruit trees exposed to different irrigation conditions. Based on the analysis results of DANS index, the distribution of it was fairly similar to the distribution of the canopy temperature and inversely proportional to the amount of supplied water. In conclusion, it was possible to estimate the water stress conditions of target crops by using the DANS index through the GPR calibration and aerial thermal image processing methods which are proposed in this study. As such, the UAV based crop water stress estimation methods were continuously studied by considering aerial thermal images which can provide temperature information of target crops. To date, many researches are being performed to realize enhanced aerial thermal image processing techniques and develop improved calibration methods for aerial thermal images through collecting more data from fields or modifying several models. We are also progressing a study to find more robust calibration model by including both diverse aerial thermal images and environmental information. Now, we expect that our suggested methods including aerial thermal image calibration method, aerial thermal image processing techniques, and DANS water stress index have good potential for use in the estimation of water stress for fruit trees as well as other scientific fields such as nondestructive evaluation of several facilities including bridges, roads, and buildings.

Data Availability
Aerial thermal images of crops were collected at the experimental farm located in National Institute of Horticultural and Herbal Science (Rural Development Administration, Korea). Aerial thermal image processing methods were realized to precisely extract the canopy temperature of fruit trees from aerial thermal images obtained by the UAV. The subjects of the experiment were 21 apple trees, located within an apple orchard (length: 40 m, width: 3 m).

Conflicts of Interest
The authors declare that they have no competing interests.