1. Introduction
Land surface temperature (LST) is the primary driving force of the exchange between turbulent heat flux and long-wave infrared (LWIR) radiation (8–15
m) at the interface of the Earth’s surface and atmosphere. Therefore, the LST is a key parameter for the physical description of the surface energy and water balance processes at the local to global scale [
1,
2]. The potential of this crucial parameter has been repeatedly demonstrated in various thermal infrared-based studies and applications, such as evapotranspiration [
3,
4], hydrological modelling [
5], vegetation monitoring [
6], ‘urban heat island and urban development’ [
7,
8,
9], climate change and weather conditions [
1,
10], agriculture [
3,
11], and the monitoring of land use changes in wetlands [
12]. In agricultural applications, precision farming has increasingly required thermal remote sensing techniques to detect water-stressed crops [
3,
13,
14], plant diseases [
13,
15], and for irrigation management [
13,
14].
However, LST significantly varies, both temporally and spatially (about 10 K), during the daytime depending on the present weather conditions and the solar irradiation that warms the land surface [
1]. Surfaces with high spectral differentiated emissivities, such as bare soil surfaces and rocks (inhomogeneous heating) and sealed areas (‘urban heat islands’) have shown particularly large surface temperature amplitudes in the diurnal cycle [
1,
16,
17]. Water-stressed vegetation also depicts this phenomenon, due to the rapidly changing transpiration rate resulting in a strongly fluctuating plant canopy temperature [
18]. Therefore, high repetition rates (depending on water supply and weather conditions) or even continuous temperature measurements are required to account for the detailed and accurate thermal mapping of water-stressed crops or in irrigation management.
Satellite-based LST retrieval methods have been reviewed and compared in several case studies [
1,
19,
20], but depending on the retrieval algorithm, large discrepancies between the retrieved LST and ground truth data can occur. Even the same sensor, such as a high precision satellite instrument (i.e., the Spinning Enhanced Visible and Infrared Imager—SEVIRI—aboard Meteosat-9) using different LST algorithms has featured discrepancies of about 6 K [
21]. The reason for this is often an inaccurate or even erroneously estimated surface emissivity [
1], which has a great impact on the accuracy of LST [
19,
20]. The lower the emissivity, the higher the impact of environmental thermal emissions on the LWIR radiation (
Section 2.1), which needs to be considered to avoid errors in LST retrieval (
Section 2.1.3). Moreover, the results also depend on properly performing atmospheric corrections because of changing weather conditions: A rise in relative humidity leads to decreasing transmittance, which attenuates the thermal signal’s path through the atmosphere [
1]. In terms of the great land surface heterogeneity, the central issue is an incorrectly estimated emissivity. Land surface emissivity (LSE) is highly dynamic, particularly in the mid and high latitudes with seasonally changing vegetation cover of miscellaneous landscape types, such as mixed forests, grasslands, wetlands, and in agricultural areas with short-term land use and land cover (LULC) changes [
1,
10,
12,
22]. This makes the precise retrieval of LST difficult, as LSE depends on changes in both natural and man-made land surface properties. Hence, an accurate and relatively timely LSE estimation is necessary for the retrieval of correct LST results.
In addition, a scale effect, due to different pixel resolutions when using the same LST algorithm for observations with different spatial scales at different altitudes, must be taken into account. For instance, a pixel might only contain a part of the plant at the leaf level, whereas the same pixel scene might represent a mixture of multiple information of bare soil and overlapping leaves at the canopy level [
1,
10,
17]. Resolution cells of mixed pixels increase the uncertainty if they are not taken into account. With respect to airborne and satellite sensed data, this influencing factor becomes more important, due to the lower spatial resolution at regional scales [
1,
17]. Therefore, high-contrast and highly spatially resolved thermal infrared mapping is urgently needed (e.g., by airborne campaigns) to investigate the mentioned scale effects and to overcome the large discrepancies observed between remotely sensed data and field measurements.
The use of thermal and multispectral optical sensors on board an unmanned aerial vehicle (UAV) has started with very simple configurations (e.g., [
23]). Although some parameters in UAV-based LST retrievals are often simplified due to low flight altitudes (e.g., atmospheric transmittance). Ref. [
24] performed a land use classification and assigned specific emissivity values from reviewed literature to each class in order to derive the actual temperature of different target surfaces. Sagan et al. [
25] used environmental parameters measured at nearby ground stations to correct air temperature, humidity, and emissivity using linear calibration equations. Similarly, Si et al. [
26] determined a single emissivity value that was applied to the entire study area. They measured the emissivity at a certain temperature and integrated it with the spectral response function of the thermal imager [
26]. Malbéteau et al. [
27] estimated emissivity in terms of vegetation cover, expressed by vegetation index values [
17,
19]. Some studies, however, do not take emissivity into account and therefore provide results that should be viewed as qualitative rather than quantitative [
14,
28]. To overcome the challenges mentioned, we have introduced a dual sensor system that was combined with a thermal imager that had a much higher sensor resolution than those available off-the-shelf. A multispectral optical sensor was used to provide vegetation indices to estimate emissivities, whereas, thermal images were captured almost without distortion due to the slow flight speed of the fixed-wing glider.
For LWIR applications, universal emissivity (
) for dense and healthy vegetation surfaces of
[
29] has typically been assumed, although the emissivities of different landscape types, such as forests, grasslands, and agricultural areas, are diverse [
17,
30,
31]. Additionally, the vegetation density has a direct impact on LSE; for instance, when there is a soil background within a resolution cell. Hence, this assumption is valid only under certain conditions. The emissivity also depends on the state of the plant, and can significantly decrease (e.g., for dry vegetation) [
31]. Generally, an error in
of 1% corresponds to a value of 0.75 K [
32].
In the thermal spectrum, fully vegetated canopies show L
ambertian behaviour, where angular dependence on emissivity is negligible [
33]. Accordingly, the thermal anisotropy can be disregarded for high emissivities [
33,
34]. However, the prerequisite for this is that the radiative transfer model has been applied (
Section 2.1). Nevertheless, thermal measurements are affected by water vapour absorption. Therefore, an atmospheric correction is reasonable to consider the water vapour. In particular, at higher relative humidities (>50%), it is expedient to determine tropospheric transmittance, even for a flight altitude (sensor-target distance) of 50 meters [
35]. Thus, we performed an atmospheric correction in this study.
In summary, if an incorrect specific emissivity is assumed and environmental thermal emissions (background temperature) and tropospheric transmittance are disregarded, a potential error of about 9 K can arise (
Section 3). Eventually, retrieval of the precise LST (or the kinetic temperature) is necessary; for instance, to assess evapotranspiration or to calculate agricultural indices such as the Crop Water Stress Index (CWSI). This is an indicator for the water supply status of plants, which is required in precision farming (e.g., for irrigation management) [
36]. An accurate estimation of LSE to generate precise LST maps remains a key prerequisite.
Thermal remote sensing provides complete spatially averaged and highly temporally resolved information about specific radiation characteristics to describe the intricacies of the LST, rather than point values by field measurements over large areas [
1]. Hence, we have investigated various algorithms using an LWIR thermal camera, as well as a multispectral sensor covering the visible and near infrared (VNIR) domain, to identify the method with the highest potential in terms of retrieving LST. Multispectral VNIR data were used to derive the Normalised Difference Vegetation Index (NDVI) [
37]. We chose predefined NDVI thresholds, where we constrained the emissivity per pixel to distinguish between bare soil surfaces, dense plant canopies, and mixed surface types (
Section 2.1.2).
The present research is focused on agricultural areas, which entail highly dynamic LSE values, due to the varying phenological statuses of the crops. Besides the natural ripening stages, this is mainly caused by external impacts, such as water stress, plant diseases, or frost damage. The objective of this study is to retrieve precise LST, taking into account the highly variable LSE derived by using an NDVI threshold approach, as well as the atmospheric impacts. For this purpose, we equipped a UAV with a multispectral sensor and a thermal camera to simultaneously capture VNIR and LWIR imagery of the entire study site. Based on the acquired VNIR and LWIR images, we generated corresponding ortho-mosaics to retrieve and quantify the LST values. Finally, we have developed a standardised LST retrieval algorithm which is optimised for fixed-wing UAV applications, in order to enable semi-supervised thermal mapping in agricultural areas. The study site of the UAV campaigns (April 2018) was an agricultural research field in Rheinbach (North Rhine-Westphalia), which is affiliated with the Agricultural Faculty of Bonn University in Germany.
2. Methodology of the LST Retrieval Algorithm Using UAV-Based Multispectral Data
The following sections describe the UAV-based LST retrieval methodology in detail. The various processing steps required to retrieve precise LST values are illustrated in
Figure 1:
For this approach, we developed a dual sensor system to simultaneously acquire multispectral VNIR and LWIR thermal imagery from a flight altitude of 77 meters, resulting in a spatial resolution of 7 cm per pixel (px) and 10 cm/px, respectively (
Figure 1i). Different reflectance panels were captured during the flights for subsequent calibrations of the VNIR and LWIR data. The background temperature was incorporated by using a panel of crumpled aluminium foil and quantified with a pyrometer [
29], which was sensitive at similar wavelengths as the thermal camera (
Appendix A.1). Due to its high specific heat capacity, the water temperature of a nearby pond was used as a reference for validation of the thermal observations (
Figure 1ii). The advantage of a water surface reference is its quite stable and homogeneous temperature distribution during the recording time and its constant emissivity (
Section 2.2.3). VNIR and LWIR digital numbers (DN) were then converted into reflectance and brightness temperature (BT), respectively. Further preprocessing steps that we performed on the VNIR data included reflectance calibration and band sensitivity normalisation. LWIR images were removed, regarding their contrast applying a 2D F
ourier transform (
Appendix A.2); afterwards, we conducted drift compensation and non-uniformity correction (
Figure 1iii). Subsequent atmospheric correction for the thermal data was based on the radiative transfer model (
Figure 1iv). Before constructing the ortho-mosaic, additional calibration of the lens parameters to rectify the LWIR imagery was necessary (
Appendix A.2). For both data sets, the same tools (with adapted parameters) to generate the ortho-mosaics were chosen: aligning imagery, then dense cloud and mesh building (
Appendix A.2). Then, we calculated the NDVI to constrain the LSE per pixel, based on predefined NDVI thresholds. Co-registration and resampling of the VNIR and BT maps was needed for pixelwise calculation of the stacked layer (
Figure 1vi). Finally, we computed the LSE, based on the proposed methods, and retrieved the LST, considering the specific emissivity, transmittance, and background temperature (
Figure 1vii).
2.1. Lst Retrieval Methodology
Based on the radiative transfer model (RTM), the radiation path through the troposphere can be written as [
17,
29]:
where
is the at-sensor radiance,
the surface radiance attenuated by the atmosphere, and
the upwelling radiance emitted by the atmosphere.
describes the downwelling radiance emitted by the atmosphere (e.g., by water vapour re-emission), which is reflected at the surface (
) and attenuated by the atmosphere. The atmospheric transmittance (
) is dimensionless (defined in the range of 0–1) and for the radiances, the unit W m
sr
is used. As
depends on both the wavelength and sensor-specific properties, the inverse P
lanck function must be used to convert from the raw at-sensor radiance (
) into the at-sensor brightness temperature (
). Planck’s law (1900) describes the spectral density of electromagnetic radiation emitted by a blackbody at a given temperature [
1]. However, the LWIR detector does not cover the entire thermal wavelength range. Thus, the camera manufacturer has used a modified P
lanck function [
38], which is based on the sensor-specific spectral response curve and laboratory-derived P
lanck constants (
Appendix A.1).
is automatically calculated by the thermal camera software, with parameters of transmittance and emissivity set to one.
Due to the reflections of natural surfaces (
) in the thermal domain, the radiometric temperature or brightness temperature sensed by the sensor is always lower than the equivalent measured absolute (kinetic) temperature (T). Radiant and kinetic temperatures are two different quantitative terms in thermodynamics [
34,
39,
40], which can be equal in no way, and simply provide the best approximations to describe the temperatures [
41]. The thermodynamic relationship for
can be described, by the S
tefan–B
oltzmann law (1884), as [
34,
40,
42]:
where
M is the irradiance (W m
),
is the emissivity (ranging from 0–1; dimensionless),
is the S
tefan–B
oltzmann constant (
W m
K
), and T is the absolute temperature (K).
According to [
29,
43], based on RTM (Equation (
1)) and using the S
tefan–B
oltzmann law (Equation (
2)), the LST retrieval (Equation (
3)) can be written as (
Appendix A.3):
where
(
) is the retrieved surface temperature,
the at-sensor brightness temperature,
the background temperature, and
the air temperature (K, respectively); and emissivity (
) and transmittance (
) range from 0–1 (dimensionless). From Equation (Equation (
3)), it is obvious that an atmospheric correction (
Section 2.1.1) considering
,
,
, and the corresponding emissivity are needed to retrieve a precise and correct absolute temperature (T) at the surface. At altitudes above 100 meters, we recommend considering the (dry) adiabatic lapse rate to correct the air temperature to the appropriate UAV height. The background temperature was incorporated using a panel of crumpled aluminium foil (
Appendix A.1). Equation (
3) has been described by [
43] and was used in an analogous way in [
29], where we included the parameter
to make the use of an additional commercial software for atmospheric correction unnecessary. The standard software (e.g., FLIR Thermal Studio Pro, FLIR Systems, Inc., Wilsonville, OR, USA; IRBIS 3 professional, InfraTec GmbH, Infrarotsensorik und Messtechnik, Dresden, Germay) uses Equations (
4) and (
5) (
Section 2.1.1) to perform atmospheric correction. To validate the retrieved LST (
Section 3), we used the water temperature (as measured by thermocouples) of a nearby pond as relative constant temperature reference (
Section 2.2.3). In order to account for small spatial temperature variation, we measured the air temperature at two different locations. Finally, we used the estimated emissivity (
) to retrieve the precise LST (
Section 2.1.2), as proposed by [
17,
44].
2.1.1. Impacts of Atmospheric Conditions on Thermal Measurements
In this subsection, we address several complex influencing factors that can affect thermal measurements during data acquisition: Infrared thermography must be executed in thermal equilibrium to obtain physically correct data. However, changing the ambient conditions of the camera (
FLIR Tau 2 640), can affect the thermal measurements and reduce the radiometric accuracy by up to ±5
C [
45], because uncooled thermal sensors (without thermo-electric cooler) cannot maintain a stabilised constant temperature of the microbolometer (focal-plane array, FPA). The accuracy of the measured
is dependent on the internal FPA (detector) temperature and housing temperature, as well as external effects caused by the UAV system itself. These parameters fluctuate throughout a flight (e.g., from shaded to unshaded conditions) and should, therefore, be taken into account [
46]. The issue of fluctuating detector and housing temperatures has been solved by the ‘Advanced Radiometry’ tool (FLIR) using a built-in shutter to calibrate itself. Before and after changing temperature events, flat field correction (FFC) is performed by submitting a uniform temperature (a flat field) to the detector. During the FFC process, the internal shutter of the camera is closed and quickly reopened to estimate the shutter temperature. Therefore, each FFC event was detected during the data acquisition and tagged in the recordings. Then, the
ThermoViewer 2.1.6 software (TeAx Technology GmbH, Wilnsdorf, Germany) was used to perform ‘Drift compensation’ (for temperature fluctuations) with appropriate correction coefficients. Drift of the housing temperature was also compensated by permanent cooling with an adjacently mounted fan. However, sensory effects can also lead to non-uniformity in thermal imagery, such as ‘cold corners’ (vignetting). The same software, hence, was used for a motion-based ‘Non-Uniformity Correction (NUC)’ filter, which distinguished constant errors by shifts within the frames.
Furthermore, thermal measurements are affected by water vapour absorption. Therefore, atmospheric correction in the thermal spectrum was performed to consider tropospheric conditions (
Section 1). The required parameters were obtained from the in situ weather station and a microclimatic station in the field (
Section 2.2.3). Initially, the present water vapour content (WVC) was computed using parameters such as the air temperature and relative humidity:
where
is the water vapour content (mm),
is the relative humidity (ranging from 0–1; dimensionless),
is the air temperature (
C), and
,
,
, and
are specific parameters for a temperature range of −40 to +120
C in the LWIR domain (
= 6.8455 × 10
,
= −2.7816 × 10
,
= 6.939 × 10
, and
= 1.5587) [
43,
47]. Subsequently, the transmittance was estimated using parameters such as the flight altitude (distance) and WVC:
where
is the transmittance (0–1; dimensionless), d is the distance (m), along with specific parameters such as the scaling factor of atmospheric damping (
= 1.9), the attenuation of the atmosphere without water vapour (
= 0.0066,
= 0.0126), and the attenuation of water vapour (
= −0.0023,
= −0.0067) [
43,
47]. Atmospheric damping is caused by absorption by gaseous components and the turbidity of the troposphere, which decrease the amplitude of the electromagnetic radiation [
48]. From this equation, it is obvious that the transmittance depends mainly on the distance and the WVC.
2.1.2. Lse Estimation Based on NDVI Thresholds
Remotely sensed spectral vegetation indices, such as the NDVI, allow a more representative spatial estimation of the emissivity than point measurements [
30,
37]:
where
and
R are the reflectances in the near-infrared and red spectra, respectively (−1...+1). Moreover,
thresholds provide a reasonable way to differentiate the emissivities of various land surface types. This is particularly necessary for dynamic land surfaces, such as seasonally varying natural vegetation and agricultural areas.
The first predictions of LSE from NDVI values were made in [
30,
44]. For estimating the emissivity determined by NDVI thresholds, the vegetation cover method was initially proposed in [
49,
50]. Using information of the fractional vegetation cover or vegetation proportion (
), along with knowledge of the emissivities of full vegetation and bare soil surfaces, this method produces effective LSE. The NDVI threshold method (NDVI-THM), using certain NDVI thresholds, was first introduced in [
51], where it was applied to Advanced Very High Resolution Radiometer (AVHRR) satellite data. For applications in agriculture, we identified this method as the most suitable, as it classifies the emissivity into three different surface types: bare soil (
<
), full vegetation (
>
), and mixed areas (composed of both attributes). Using information from
and NDVI histograms to identify the convenient NDVI thresholds is acutely beneficial: No accurate atmospheric correction is essential, because corrected and uncorrected NDVI are similar [
17,
52]. For known emissivities of full vegetation and bare soil surfaces,
(Equation (
7)) can be obtained, according to [
17,
52], by:
where
and
are the NDVI of bare soil and full vegetation, respectively. The
values can, then, be used to compute the emissivity of mixed pixels (
):
where
(=0.988) and
(=0.935) are the mean emissivities (in the LWIR spectrum) of dense plant canopy [
30,
31] and bare soil (of loamy to silty texture) [
31,
53],
C is the term of cavity effect due to surface roughness (
), and
is a revised parameter (d
≈ 0.01). For flat areas (
), the formula (Equation (
8)) simplifies to:
Finally, we constrain the emissivity per pixel based on certain NDVI thresholds to distinguish between three different surface types, as proposed by [
17]. The following NDVI thresholds were used to reclassify the LSE (
), based on the emissivities of bare soil, full vegetation, and mixed states:
bare soil: | | ⇒ NDVI < 0.157 | → = 0.935 |
dense canopy: | | ⇒ NDVI > 0.905 | → = 0.988 |
mixed pixels: | | | → (Equation (9)), |
where
is maintained as the global emissivity of dense canopies [
30,
31] and
is assumed for loamy to silty (moisture) soil [
31,
53]. For loamy sand soil, we suggest a lower emissivity of 0.914, as proposed by [
30]. A
value of 0.157 (standard deviation, SD = 0.227) was derived by the means of four samples from the rapeseed field studied. A
value of 0.905 (SD = 0.111) was established by the means of four rape plots (healthy varieties). However, for different crops (i.e., corn, wheat, barley, and rape), we recommend a global
value of 0.814 (as observed in summer 2018); or around 0.8, as proposed by [
17]. Nevertheless, the NDVI thresholds for soil and vegetation need to be recalculated before using NDVI-THM to estimate the LSE.
For open canopies, ref. [
30] have postulated a logarithmic relationship (
) between the mean emissivity and the mean NDVI, which fits best to NDVI values between 0.157 and 0.727 (Equation (
10)):
with
a = 1.0010 and
b = 0.047 (level of significance: 0.99). To analyse the results of the emissivity approach from information derived from
(Equation (
9)), we used the logarithmic relationship (Equation (
10)). Finally, we performed a comparison of both approaches to determine their selectivities, as well as to study their benefits for applications in agriculture (
Section 3).
2.1.3. Consequences of Omitting Essential Parameters
Subsequent assumptions were made to emphasise the importance of considering essential parameters, such as transmittance, air temperature, emissivity, and background temperature. In our first assumption, we set the emissivity to one (
1), which simplifies the formula (Equation (
3)) to:
In our second assumption, which presupposes a completely transmissive atmosphere, we set the transmittance to one (
1). This simplifies the same initial formula to:
Based on these limiting assumptions (Equations (
11) and (
12)), we made the following two case distinctions: Case i: LST(
) includes the atmospheric transmittance, but does not take into account the emissivity or background temperature (Equation (
11)); and case ii: LST(
) includes the emissivity and background temperature, but omits the transmittance (Equation (
12)). In summary, we have highlighted the potential implications of omitting these parameters, which should be considered in order to minimise the uncertainties in retrieving precise LST values. The deviations between the uncorrected
and LST(
) or LST(
) are shown in the results (
Section 3).
2.2. Study Site and Data Sets
2.2.1. Agricultural Research Campus
The study site of the UAV campaigns was the northern part of agricultural research fields of ‘Campus Klein-Altendorf’ (CKA) in Rheinbach (North Rhine-Westphalia), which is affiliated with the Agricultural Faculty of Bonn University, Germany (
Figure 2). The total agricultural area of CKA was about 181 ha, mainly cultivated with barley, wheat, corn, and sugar beet; moreover, rapeseed, potatoes and renewable resources (e.g., Miscanthus) were also cultivated. CKA had two weather stations that continuously recorded meteorological data, such as air temperature and humidity, as well as the surface temperature in the northern and southern part of the campus. In the fields, air temperature and humidity were measured continuously by data loggers (
Section 2.2.3). This study was focused on the rapeseed field within the ‘Common Experiment’ (
Figure 2) during the springtime in 2018. The maintenance of crop rotation and working on homogeneous fields was part of this experiment. The rapeseed field experiment consisted of 63 plots (and border plots), each with dimensions of 3 × 6 meters and planted with 10 different genotypes with a minimal repetition of 6 plots (the present genotypes were Alaska, Pirola, DH5, Markus, Expert, Olympiade, Wotan, Major, Groß Lüsewitzer, and Sensation NZ).
2.2.2. Fixed-Wing Uav System
A fixed-wing UAV (
Figure 3) was used to map the heterogeneous and diverse agricultural cropping systems and to efficiently cover the large fields of the campus.
The special features of the UAV are the novel sensor platform (
Figure 3) and a slow flight speed over ground (
Appendix A.4), which is particularly important for thermal mapping, in order to avoid blurry images. The UAV was mounted with a dual sensor system, consisting of both a multispectral VNIR sensor (Parrot Sequoia; Parrot Drones SAS, Paris, France) and a LWIR thermal camera (
FLIR Tau 2 640), attached to a
ThermalCapture 2.0 module (TeAx Technology) to simultaneously capture the images (
Table 1). The accuracy of the thermal camera (
) under stable ambient conditions was
= ±0.50
C, which has been laboratory-confirmed by [
45]. However, in changing ambient conditions, the radiometric accuracy can be reduced by up to ±5.0
C [
45].
An irradiance-sensing element, belonging to the VNIR sensor, was used to account for the present lighting conditions. This sensor was fixed on the UAV and oriented towards the sky. Both sensors covered the optical spectrum, separated into four bands (centre wavelength): green/G (550 nm), red/R (660 nm), red edge/RE (735 nm), and NIR (790 nm). Moreover, the sensor platform was equipped with an internal fan to cool the devices and to provide a stable ambient temperature, which is relevant for proper thermal imaging. An advanced autopilot (Pixhawk 4) connected to a GPS antenna for positioning and a pitot tube (a device measuring the flight speed with respect to wind speed) on the wing were built-in to automatically fly the programmed routes. Generally, we conducted the UAV campaigns once or twice times per month during the growing season. The flights were split into three parts (approximately 60 ha per section), due to the limited power supply of the UAV batteries, which were generally scheduled at 10:00 (focused study field), 12:00, and 14:00 (CET), to minimise the influence of shadows in the imagery. Every flight took approximately 30–45 minutes.
2.2.3. Meteorological Data and Reference Measurements
A weather station (
GWU-Umwelttechnik GmbH) was situated in the northern part of the campus (near the pond), which recorded the following parameters (at 2 m above the ground): wind speed, air temperature (
), relative humidity (
), surface temperature (
) at 0.05 m, and net radiation (
Table 2 and
Table 3, respectively). Within the rapeseed field, a microclimatic station continuously recorded (at 1 min intervals) the air temperature and relative humidity (
HMP110) as shown in
Table 3. A pull-up resister (
DS18B20) measured the (bare soil) surface temperature with an accuracy of ±0.5
C. In addition to the water surface temperature of the pond,
was used to validate the retrieved LST from the thermal camera (
Section 3).
As the weather conditions were quite stable (
Table 2) and there were no noticeable fluctuations in
or abruptly changing plant states (such as crop water stress), no further reference measurements (aside from the water temperature) were carried out with thermocouples to validate the surface temperatures [
29]. However, weather conditions have a direct impact on the surface temperature: increasing
will directly increase
in analogical degree, but also change the net radiation, affecting the surface temperature. The interdependence of
and
as an approach to surface temperature correction (
) has been proposed by [
29]:
where
is the corrected surface temperature,
,
= <
>, and
is the mean of <
> (
C). The correction of
compensates for the noise caused by drifts in air temperature (
Figure 4).
This approach (Equation (
13)) was used to account for small discrepancies (
) in the thermal measurements caused by the impact of <
>. We calculated the differences of
and
, which were defined as
. There were no significant changes in <
> and
during the campaign and, so,
= 0.28 ± 0.1
C was small.
Due to stable weather conditions (
Table 2), the averaged values of
(from before, during, and after the flights) were used (
Appendix A.1) to consider the background temperature [
29]. Moreover, for cloudy conditions (diffuse scattering), the background temperature is quite similar to the measured air temperature (
Table 3). Flight 3 showed a much lower background temperature (
), due to (‘cold’) clear sky conditions. The transmittance (
) was always about 0.95. Both
and
were used to perform atmospheric correction.
In addition, during the flights, data at a nearby pond were captured, which were used as reference for the temperature validation (
Section 3). The Normalised Difference Water Index (NDWI) was applied to locate the surface of the pond [
54,
55]:
where
G and
are the reflectances in the green and near-infrared spectra, respectively (−1...+1). According to [
56], we used an
threshold of
≥ 0.3 to detect the water surface and to reclassify the corresponding emissivity
(
= 0.985) [
31,
34]. Water surfaces almost totally absorb LWIR at a layer thickness of >1 mm and show a constant emissivity [
40].
In terms of the water surface temperature (
) measurements, we constructed a floating triangle out of bamboo poles (with a side length of 3 m each) for the thermocouples, which were mounted at each corner and in the middle. The thermocouples were placed at depths of 3.0–30 mm (the dimensions of a thermocouple element is 5.1 × 33 mm).
was measured directly under the surface layer by four thermocouples, using a data logger with an accuracy of ±0.36
C. The mean temperature of the four thermocouples was used to validate the retrieved camera temperature (
). According to Equation (
3), the captured thermal images of the pond were used to retrieve
. Before and after the flight over the Common Experiment (
Figure 2), each set of 24 suitable thermal images of the pond was selected for validation. Both
and
were used to calculate the temperature deviations (
) during the flight.
3. Results
The CIR composite (
Figure 5a) enhanced the visibility of different varieties for the rape field studied:
Similar reddish colours showed the same varieties, depending on their biomass, and dark green plots were plants damaged by frost due to cold spells in the springtime. These plants belonged to the same genotype and turned out to be less resilient against frost. Plants affected by frost were inhibited in their growth or entirely perished. Striking greenish patches depicted bare soil areas around the plots. These patterns continued in the subsequent maps. The NDVI map (
Figure 5b) was generated to find certain NDVI thresholds of bare soil (
= 0.157; SD = 0.227) and full vegetation (
= 0.905; SD = 0.111), which were used to compute the vegetation proportion (
) and to estimate the LSE (
Section 2.1.2). The
map (
Figure 5c) was used to compute the LSE (Equation (
9)), whereas LSE’ was calculated based on the logarithmic relationship (Equation (
10)). High NDVI values represent a healthy and dense canopy; corresponding to
= 1 for full vegetation. Both the NDVI and the
maps confirm the lack of biomass for the withered varieties, resulting in low NDVI. The latter map even indicated a marginal vegetation proportion (
= 0) for the affected plots, similar to that of bare soil surfaces.
The LSE (
Figure 6a) and LSE’ (
Figure 6b) maps were computed using certain NDVI thresholds for bare soil, dense canopy, and mixed states (
Section 2.1.2). For mixed pixels, LSE was estimated by information from
(Equation (
9)) and LSE’ by using the logarithmic relationship (Equation (
10)). Comparison of LSE and LSE’ was performed to study the selectivity of both approaches and to evaluate the accuracy of the NDVI-THM approach (
Figure 6c). The minimum emissivity was 0.935, representing loamy to silty soil, and the maximum emissivity was 0.988, representing dense canopy. The LSE (
Figure 6a) map, derived by the
approach, shows detailed information of the open canopy, such as withered plants and plot edges (white areas with an emissivity of about 0.96). The LSE’ map (
Figure 6b), based on the logarithmic relationship (
’) between the mean emissivity and the mean NDVI, fits for NDVI values between 0.157 and 0.727 better than for NDVI thresholds of up to 0.905 (i.e., for dense-canopied rapeseed).
In contrast, the
approach could be adjusted to determine the certain NDVI values of different cereals. The comparison (
Figure 6c) of both approaches shows a selectivity of dLSE = 0.03, especially for plants affected by frost and plot edges (red areas). This fact could be an indication of a discrepancy in the estimation of lower emissivities (i.e., for dry vegetation). For dense canopy and bare soil, dLSE asymptotically approached zero, which confirms the accuracy of the estimated LSE. However, we have stated that the NDVI-THM approach cannot be used to estimate the emissivity of senescent or withered plants, as the NDVI of dry vegetation (which is considerably affected by water stress) has a similar emissivity spectrum as that of bare soil [
57,
58]. The emissivity of dry vegetation can reach very low values, compared to healthy plants, which leads to LST values that are considerably underestimated; without additional consideration of meteorological parameters, a potential error of about 9 K can occur (
Section 4).
The
map (
Figure 7a) clearly distinguishes colder areas of the dense canopy and warmer spots (up to 4.8 K) of the rape varieties affected by frost, as well as bare soil surfaces. As the measured at-sensor brightness temperature does not take into account the surface emissivity and current atmospheric conditions,
should only be used for surfaces with high emissivities and relative temperature observations under similar weather conditions. The LST retrieval (
Figure 7b) included the estimated LSE (
Figure 6a), as well as the present atmospheric conditions. The minimum temperature of the dense canopy was 13.5
C and the maximum temperature of the bare soil and dried-up plants was 19.2
C. The difference between LST and
(
Figure 7c) showed a maximum deviation (dLST) of 1.0 K for bare soil surfaces and varieties that had suffered from frost. The accuracy of the LST retrieval is strongly dependent on the specific surface emissivity and differs substantially with decreasing emissivity: Without consideration of mixed-resolution cells from soil and vegetation, the LST will be underestimated for affected pixels by at least 1 K (
Figure 7c). Hence, we recommend using an NDVI thresholding method (
Section 2.1.2) for the emissivity estimation, in order to consider mixed pixels and to minimise the uncertainties.
The final LST retrieval (
Figure 7b) depicts almost the same temperature distribution as the
map: warmer spots (up to 5.7 K) showed frost-damaged plants and bare soil, and dense canopy occurred as colder areas. Compared to the CIR map, with specific reddish colours for same varieties, the temperature variation within the plots was larger than for different varieties. Therefore, the varieties were not distinguishable by certain LST thresholds.
Moreover, we had nearly constant weather conditions during the campaigns (
Section 2.2.3) and no noticeable shifts in the plant state. Thus, the temperature variance of the varieties was generally caused by different evapotranspiration, plant water content, and soil water supply. In our case, the differing canopy temperature within the plots was mainly related to the diverse emissivities of withered and healthy plants.
Finally, we analysed the impacts of the meteorological parameters on thermal measurements. Different assumptions were made to emphasise the effects of when essential parameters, such as transmittance and environmental thermal emissions, were not taken into account in LST retrieval (
Section 2.1.3). Each case was compared to the uncorrected
(
Figure 7a) and the modified LST retrievals (Equations (
11) and (
12)) to depict the consequences of inaccurately retrieved LSTs: Case i: Atmospheric transmittance included in LST(
), but without consideration of the background temperature and specific emissivity (Equation (
11)); and case ii: The background temperature and specific emissivity included in LST(
), but without consideration of transmittance (Equation (
12)). For
without consideration of the atmospheric transmittance, a temperature deviance of up to 0.4 K occurred (
Figure 8a). This was mainly caused by water vapour absorption in the thermal domain. In contrast, omitting the background temperature and specific emissivity led to a deviance of up to 0.6 K (
Figure 8b). This result was derived during cloudy sky conditions (
C). In a modified scenario, with respect to the background temperature in clear sky conditions, the temperature deviation was up to 2.9 K (
Figure 8c). In all scenarios, the largest deviation was found for bare soil areas. During clear sky conditions,
had minimal temperature −40
C and maximal temperature 10
C, in the case of a diffusely scattering sky [
32]. Diffuse scattering (in overcast conditions) leads to almost identical measurements between the air temperature and the observed background temperature (
). However, the background temperature should be taken into account to minimise uncertainties and to prevent the LST from being underestimated. In conclusion, environmental thermal emission (if
< 1) is the most influential meteorological parameter in thermal measurements.
The NDWI map (
Figure 9a) was generated to detect the pond’s surface and to reclassify the corresponding emissivity (
) of the water using a determined NDWI threshold (
Section 2.2.3). Based on the emissivity map (
Figure 9b), the mean pond surface temperature (
) was retrieved (
Figure 9c). Before and after the flight over the agricultural experiment site (flights 1 and 3, respectively), we retrieved 24 suitable LST images (
) to validate the thermal data using the measurements with the pond thermocouples (
). Both
and
were used to calculate the mean temperature deviations (
) to assess the accuracy of retrieved LST.
The mean absolute error (MAE) was used to compute the errors of the temperature comparisons [
59]. Flight 1 (
Figure 10a) had an error of
= 0.53 K, with a coefficient of determination of
= 0.985 (
n = 24). For flight 3 (
Figure 10b), the MAE was
= 0.54 K and
= 0.964 (
n = 24).
The results of the performed validation, especially for flight 1, confirmed the high significance of the performance of the retrieval algorithm (Equation (
3)). During flight 3, the correlation was somewhat less significant in comparison to flight 1.
The temperature validation of the rapeseed experiment (
) also confirmed the high performance of the LST retrieval (
Figure 11). The MAE was
= 0.41 K and
= 0.879 (
n = 27).
For the UAV campaigns using the proposed LST retrieval algorithm (Equation (
3)), we found the accuracy of the thermal measurements to be around MAE = 0.5 K.
4. Discussion
In this study, we have developed an UAV-based LST retrieval algorithm (Equation (
3)) suitable for conducting thermal campaigns, especially in agricultural areas. Agricultural applications increasingly demand precise and current LST maps for irrigation scheduling, evapotranspiration and drought stress monitoring, and plant disease detection. A true physical temperature, as calculated by the LST retrieval of crops, is essential for calculating indices such as the Vegetation Health Index (VHI) or CWSI [
13]. The accuracy of retrieved LST depends mainly on properly estimating the emissivity of the surface type (LSE). Thus, we have used the proposed NDVI threshold method (
Section 2.1.2) to distinguish the emissivities in three classes: bare soil (
), dense canopy (
), and mixed states (Equation (
9)). For dense canopy, we suggest a global
of 0.814 (as observed in summer 2018) or around 0.8 (as proposed by [
17]), which was averaged for different crops (i.e., corn, wheat, barley, and rape). For bare soil, we derived
to be 0.157, in agreement with [
30], or around 0.2 (as proposed by [
17]). In the case of higher sand content, lower NDVI values arise.
The determined emissivity of dense canopy (
= 0.988) agreed with frequently cited values, as in [
27,
29,
31,
60]. However, the plants must be healthy, as the emissivity decreases significantly in the case of withered crops. In the study site, soil surfaces predominantly consisted of loamy to silty textures. Thus, we maintained an emissivity of bare soil of
= 0.935 (as proposed by [
31,
53]). For a sandy soil texture, we suggest a lower emissivity, such as 0.914 (as proposed by [
30]). Regardless of the emissivity of mixed resolution cells, the LST was underestimated for the affected pixels by at least 1 K (
Figure 7c). The lower the emissivity, the higher the impact of the background temperature (
) on the retrieved LST. Hence, we recommend atmospheric correction, especially in (‘cold’) clear sky conditions, in which observation periods commonly occur to ensure the best radiometric image quality. However, this method reaches its limits in the case of lower emissivities (e.g.,
), such as those for dry vegetation [
31,
34]. Disregarding essential meteorological parameters (
Section 2.1.3) and assuming an overall emissivity for dense canopy (for instance, of
= 0.985), a potential error of about 9 K occurs.
For data validation, we continuously measured the water temperature of a nearby pond by thermocouples, which we compared with each set of 24 suitable retrieved LST images of the pond from the flights before and after observing the experiment. Water bodies provide a convenient opportunity for reference measurements, due to their relatively homogeneous and constant surface temperature. For the pond, we assumed an emissivity of
= 0.985 (according to [
31,
34]). The mean emissivity (in the LWIR range of the camera) of distilled water (
= 0.986), obtained from the emissivity library of the Moderate Resolution Imaging Spectrometer (MODIS) confirms this assumption [
61].
The advantages and disadvantages of conducting thermography under homogeneously cloudy or clear sky conditions remain to be discussed: Clear sky conditions provide the most contrasting results; and cloudless weather is indispensable for the validation of airborne or satellite data with UAV observations. However, the lower the emissivity, the higher the reflection of the surface. Therefore, the impact of environmental thermal emissions on bare soils and dry vegetation is the greatest, which must be taken into consideration. Cloudy conditions have a strong impact on the quality of the captured thermal imagery, since they reduce the contrast and radiometric resolution [
13]. However, the impact of environmental thermal emissions is low (isotropic radiation) and the background temperature is quite similar to the air temperature. Scenarios with the biggest drawbacks are rapidly changing air and surface temperatures or cloud shadows in the recordings. In order to ensure the best lighting conditions, we recommend thermal measurements in clear sky or, at least, in homogeneously cloudy conditions.
In conclusion, the validations showed MAEs of
= 0.53 K (
= 0.985) and
= 0.54 K (
= 0.964), respectively. These results confirm the highly significant performance of the retrieval algorithm (Equation (
3)). For the thermal measurements, we stated an MAE of about 0.5 K.