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

Authors: Yunhyeok Han; Barnabas Abraham Tarakey; Suk-Ju Hong; Sang-Yeon Kim; Eungchan Kim; Chang-Hyup Lee; Ghiseok Kim
Journal: Journal of Sensors
 (2021)

---

## Abstract

Crop monitoring is a very important area of precision agriculture and smart farming. Through an accurate monitoring, it is possible to more efficiently manage the irrigation, fertilization, and pest control. In this study, we propose aerial thermal image calibration method and thermal image processing techniques to analyze the water stress level of fruit trees under different irrigation conditions. The calibration was performed using Gaussian process regression, and it was demonstrated as an appropriate regression method because it satisfied all requirements including the residuals’ normality, independence, and homoscedasticity. In addition, an appropriate image processing technique was necessary to selectively extract only the canopy temperature from the aerial thermal images, while excluding irrelevant elements such as the soil and other objects. For the image processing techniques, three methods (Gaussian mixture model, Otsu binarization algorithm, and Otsu binarization algorithm after Gaussian blurring) were employed. The Gaussian mixture model provided the highest accuracy and stable results for the extraction of the canopy temperature. After the aerial thermal images were subjected to calibration and image processing, the degree above nonstressed canopy (DANS) water stress index was calculated for the fruit trees under different water supply conditions. The distribution of the DANS water stress index was similar to the distribution of the canopy temperature and inversely proportional to the amount of supplied water content. Therefore, we expect that the DANS water stress index, calculated using the calibration and image processing techniques proposed in this study, can be a reliable measure for the estimation of the water stress of crops for the application of aerial infrared techniques to remote sensing.


---

## Body

## 1. 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–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 high-spatial-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 differences 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.


## 2. 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.Figure 1
Schematic of the aerial thermal image calibration experiment.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 Figure2. (2) The blackbody’s temperature was measured on the ground, using resistance temperature detector (RTD) 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.Figure 2
Aerial thermal image of the blackbody captured at an altitude of 10 m.The regression was based on Eq. (1), which predicts the difference between the blackbody temperature measured by the RTD temperature sensors (TbG) and blackbody temperature measured by the UAV (TbA) with independent variables of altitude (Alt), atmospheric temperature (Temp), relative humidity (Humd), and solar irradiance (Irrad),
(1)TbG‐TbA=funAlt,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 data [25]. The conditional distribution of y∗ can be expressed as a multivariate normal distribution,
(2)y→∗y→,X,X∗~Nμ∗,Σ∗,(3)μ∗=KX,X∗KX,X+σ2I−1∗y→,(4)Σ∗=KX∗,X∗−KX,X∗KX,X+σ2I−1KX,X∗,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,
(5)KRQx,x,=σf21+x−x′22αl2−α,where three hyperparameters, σf2 (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 and test sets were set to 80% and 20%, respectively. For each set, we evaluatedR2 (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 the graph between fitted values and residuals. A Durbin–Watson (DW) test was performed for the analysis of independence between residuals [28]. The GPR model and its evaluation methods were developed using a commercial programming tool (MATLAB ver. R2016, MathWorks Inc., Natick, MA, USA).


### 2.2. Aerial Thermal Image Processing

As shown in Figure3, 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.Figure 3
Aerial thermal image of the fruit tree captured at an altitude of 10 m. The red box is used when the canopy temperature is manually measured.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 Figure4, 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 distributions 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.Figure 4
Histogram of the Gaussian mixture model for the aerial infrared thermal image.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,
(6)costfunction=p1×V1+p2×V2,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 changesp1 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).


### 2.3. 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 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,
(7)DANSh=Tch−TcNSh,where Tc is the leaf temperature at given time h and TcNS 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.


## 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.Figure 1
Schematic of the aerial thermal image calibration experiment.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 Figure2. (2) The blackbody’s temperature was measured on the ground, using resistance temperature detector (RTD) 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.Figure 2
Aerial thermal image of the blackbody captured at an altitude of 10 m.The regression was based on Eq. (1), which predicts the difference between the blackbody temperature measured by the RTD temperature sensors (TbG) and blackbody temperature measured by the UAV (TbA) with independent variables of altitude (Alt), atmospheric temperature (Temp), relative humidity (Humd), and solar irradiance (Irrad),
(1)TbG‐TbA=funAlt,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 data [25]. The conditional distribution of y∗ can be expressed as a multivariate normal distribution,
(2)y→∗y→,X,X∗~Nμ∗,Σ∗,(3)μ∗=KX,X∗KX,X+σ2I−1∗y→,(4)Σ∗=KX∗,X∗−KX,X∗KX,X+σ2I−1KX,X∗,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,
(5)KRQx,x,=σf21+x−x′22αl2−α,where three hyperparameters, σf2 (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 and test sets were set to 80% and 20%, respectively. For each set, we evaluatedR2 (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 the graph between fitted values and residuals. A Durbin–Watson (DW) test was performed for the analysis of independence between residuals [28]. The GPR model and its evaluation methods were developed using a commercial programming tool (MATLAB ver. R2016, MathWorks Inc., Natick, MA, USA).


## 2.2. Aerial Thermal Image Processing

As shown in Figure3, 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.Figure 3
Aerial thermal image of the fruit tree captured at an altitude of 10 m. The red box is used when the canopy temperature is manually measured.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 Figure4, 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 distributions 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.Figure 4
Histogram of the Gaussian mixture model for the aerial infrared thermal image.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,
(6)costfunction=p1×V1+p2×V2,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 changesp1 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).


## 2.3. 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 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,
(7)DANSh=Tch−TcNSh,where Tc is the leaf temperature at given time h and TcNS 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.


## 3. 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 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 Figure5, 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.Figure 5
Scatter plots between (TbG–TbA) and altitude, solar irradiance, atmospheric temperature, and relative humidity.Figure 6
Dependence of the variables. A high linearity is observed 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 Table1. 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.Table 1
VIFs for the four variables.
Altitude (m)Solar irradiance (W/m2)Atmospheric temperature (°C)Relative humidity (%)Variance inflation factor (VIF)1.03.0488.4113.885Table 2
VIFs for the three variables without atmospheric temperature.
Altitude (m)Solar irradiance (W/m2)Relative humidity (%)Variance inflation factor (VIF)1.01.1951.195


### 3.2. 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. TheR2 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 R2 for the training and test sets.Figure 7
Expected response vs. actual response: (a) training set; (b) test set.
(a)(b)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 Figure8, 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.Figure 8
Fitted values vs. residuals of the GPR for the linearity and homoscedasticity of the residuals: (a) training set; (b) test set.
(a)(b)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/m2, 45.7%, and 35.9°C, respectively.Figure 9
Corrected thermal image after GPR calibration.


### 3.3. Canopy Temperature for Image Processing Methods

Three types of aerial thermal image processing method were applied to precisely extract the canopy temperature of the fruit trees. These image processes were carried out after we cropped the tree of interest from the entire aerial thermal images, as shown in Figure10. 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.Figure 10
Aerial thermal images obtained by different processing methods: (a) without processing; (b) Gaussian mixture model; (c) Otsu binarization method; (d) Otsu binarization with Gaussian blurring.
(a)(b)(c)(d)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–14) under 50% ET, four trees (Nos. 15–18) not subjected to irrigation, and last three trees (Nos. 19–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 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.Figure11 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 quantitatively extracted the canopy temperature of the crops to make it easy and clear for the further data processing.Figure 11
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. They were taken on different four days: (i) August 7th, (ii) August 13th, (iii) August 22nd, and (iv) August 29th.


### 3.4. DANS Index for Different Irrigation Conditions

Table3 presents the canopy temperature of the target crops by day. The average canopy temperatures of the target crops on August 7th, 13th, 22th, and 29th were 46.02, 52.48, 43.62, and 32.79°C, respectively. The highest canopy temperatures were 48.3, 55.8, 46.2, and 36.5°C, while the lowest temperatures were 42.5, 50, 40.7, and 29.5°C, respectively. The standard deviations of the canopy temperature distributions were 1.57, 1.15, 1.53, and 1.92°C, respectively. The day with the highest average temperature was August 13th, while the day with the most severe temperature deviation was August 29th. The average atmospheric daytime temperatures on August 7th, 13th, 22th, and 29th were estimated to be 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.Table 3
Canopy temperature of the target crops on different days.
Temperature (°C)August 7, 2019August 13, 2019August 22, 2019August 29, 2019Average46.0252.4843.6232.79Highest48.3055.8046.2036.50Lowest42.5050.0040.7029.50Standard deviation1.571.151.531.91In Figure12, 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 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.Figure 12
Comparisons of the DANS indices of apple trees with different irrigation conditions after calibration with the GPR method: (a) August 7th, 2019; (b) August 13th, 2019; (c) August 22th, 2019; (d) August 29th, 2019.
(a)(b)(c)(d)


## 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 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 Figure5, 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.Figure 5
Scatter plots between (TbG–TbA) and altitude, solar irradiance, atmospheric temperature, and relative humidity.Figure 6
Dependence of the variables. A high linearity is observed 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 Table1. 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.Table 1
VIFs for the four variables.
Altitude (m)Solar irradiance (W/m2)Atmospheric temperature (°C)Relative humidity (%)Variance inflation factor (VIF)1.03.0488.4113.885Table 2
VIFs for the three variables without atmospheric temperature.
Altitude (m)Solar irradiance (W/m2)Relative humidity (%)Variance inflation factor (VIF)1.01.1951.195


## 3.2. 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. TheR2 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 R2 for the training and test sets.Figure 7
Expected response vs. actual response: (a) training set; (b) test set.
(a)(b)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 Figure8, 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.Figure 8
Fitted values vs. residuals of the GPR for the linearity and homoscedasticity of the residuals: (a) training set; (b) test set.
(a)(b)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/m2, 45.7%, and 35.9°C, respectively.Figure 9
Corrected thermal image after GPR calibration.


## 3.3. Canopy Temperature for Image Processing Methods

Three types of aerial thermal image processing method were applied to precisely extract the canopy temperature of the fruit trees. These image processes were carried out after we cropped the tree of interest from the entire aerial thermal images, as shown in Figure10. 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.Figure 10
Aerial thermal images obtained by different processing methods: (a) without processing; (b) Gaussian mixture model; (c) Otsu binarization method; (d) Otsu binarization with Gaussian blurring.
(a)(b)(c)(d)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–14) under 50% ET, four trees (Nos. 15–18) not subjected to irrigation, and last three trees (Nos. 19–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 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.Figure11 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 quantitatively extracted the canopy temperature of the crops to make it easy and clear for the further data processing.Figure 11
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. They were taken on different four days: (i) August 7th, (ii) August 13th, (iii) August 22nd, and (iv) August 29th.


## 3.4. DANS Index for Different Irrigation Conditions

Table3 presents the canopy temperature of the target crops by day. The average canopy temperatures of the target crops on August 7th, 13th, 22th, and 29th were 46.02, 52.48, 43.62, and 32.79°C, respectively. The highest canopy temperatures were 48.3, 55.8, 46.2, and 36.5°C, while the lowest temperatures were 42.5, 50, 40.7, and 29.5°C, respectively. The standard deviations of the canopy temperature distributions were 1.57, 1.15, 1.53, and 1.92°C, respectively. The day with the highest average temperature was August 13th, while the day with the most severe temperature deviation was August 29th. The average atmospheric daytime temperatures on August 7th, 13th, 22th, and 29th were estimated to be 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.Table 3
Canopy temperature of the target crops on different days.
Temperature (°C)August 7, 2019August 13, 2019August 22, 2019August 29, 2019Average46.0252.4843.6232.79Highest48.3055.8046.2036.50Lowest42.5050.0040.7029.50Standard deviation1.571.151.531.91In Figure12, 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 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.Figure 12
Comparisons of the DANS indices of apple trees with different irrigation conditions after calibration with the GPR method: (a) August 7th, 2019; (b) August 13th, 2019; (c) August 22th, 2019; (d) August 29th, 2019.
(a)(b)(c)(d)


## 4. 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.



---

*Source: 5537795-2021-05-11.xml*