Drift Correction of Lightweight Microbolometer Thermal Sensors On-Board Unmanned Aerial Vehicles

: The development of lightweight sensors compatible with mini unmanned aerial vehicles (UAVs) has expanded the agronomical applications of remote sensing. Of particular interest in this paper are thermal sensors based on lightweight microbolometer technology. These are mainly used to assess crop water stress with thermal images where an accuracy greater than 1 ◦ C is necessary. However, these sensors lack precise temperature control, resulting in thermal drift during image acquisition that requires correction. Currently, there are several strategies to manage thermal drift effect. However, these strategies reduce useful ﬂight time over crops due to the additional in-ﬂight calibration operations. This study presents a drift correction methodology for microbolometer sensors based on redundant information from multiple overlapping images. An empirical study was performed in an orchard of high-density hedgerow olive trees with ﬂights at different times of the day. Six mathematical drift correction models were developed and assessed to explain and correct drift effect on thermal images. Using the proposed methodology, the resulting thermally corrected orthomosaics yielded a rate of error lower than 1 ◦ C compared to those where no drift correction was applied.


Introduction
World agriculture faces three major challenges that represent an apparent contradiction: to feed a growing population, to contribute to the reduction of rural poverty, and to manage the natural resource base [1,2]. Precision agriculture is believed to be an efficient method of crop production because it is accurate, inputs are optimized leading to reduced costs and environmental impact, and because it provides an audit trail that consumers and legislation require [3]. Precision agriculture emphasizes spatial-temporal data analysis and management jointly rather than singularly [4], as well as requiring a detailed description of canopy status and its variation in the field during the growth cycle [1].
Remote sensing methods have been demonstrated to be very useful in monitoring large areas while remaining cost effective [5]. Traditional remote sensing techniques have used manned aerial or satellite platforms to measure canopy reflectance in the electromagnetic spectrum range from 400 to 2500 nm [6]. These platforms have a temporal and spatial resolution that limits their utility in

Materials and Methods
The presented study was carried out in Córdoba, Spain (north latitude 37 • 56 05 , west longitude 4 • 42 59 , WGS 84) in June 2016, using a 10 ha orchard of high-density hedgerow olive trees (Olea europaea L. cv Arbequina). Figure 1 shows the development of olive trees during the growing season. Being a typical Mediterranean region, the climate is characterized by warm, dry summers and cool, wet winters with an average annual rainfall equal to 180 mm.

UAV Campaigns
The unmanned aerial vehicle (UAV) used was an MD4-1000 multi-rotor drone (Microdrones GmbH, Siegen, Germany). This UAV is a quadcopter with an entire carbon design. The system has a maximum payload equal to 1.2 kg. It uses 4 × 250 W gearless brushless motors and reaches a cruising speed of 15.0 m/s. The UAV was equipped with a Gobi 640-GiGe thermal sensor (Xenics nv, Leueven, Belgium), which is an uncooled long-wave infrared (LWIR) thermal sensor delivering raw digital images at 16 bits of sensor calibrated radiance with a dynamic range from −20 • C to 120 • C and a spectral resolution of 0.05 • C. It has a focal length equal to 18 mm and operates in a spectral band range from 8 µm to 14 µm. Registered images have a dimension equal to 640 × 480 pixels and a pixel pitch of 17 µm. Moreover, it has an onboard image processing system to perform a non-uniformity correction, an auto offset, and an auto gain. However, the continuous changing conditions in which the sensor operates cause temperature values to degrade throughout the UAV flight although the image processing system is operating. In this manuscript, the NUC has been deactivated, obtaining a set of images without any compensation. A stabilization procedure on the thermal sensor was conducted before each UAV flight as described in Berni et al. [36]. The thermal sensor was pre-heated for twenty minutes on field before each UAV flight to stabilize its internal temperature. The sensor was connected to a stick PC Asus QM1 (Asustek Computer Inc., Taiwan, China) to store the images via an Ethernet port. Sensor weight totaled 710 g with flight duration equal to 20 min.
UAV flights were performed on 1 June 2016 at 120 m above ground level with a ground sample distance (GSD) equal to 11.3 cm (Figure 2). Side and forward lap settings were 80% and 70%, respectively, with 632 images registered. Five aluminum disks were placed on the plot as ground control points (GCPs), one in each corner and the other in the center of the study area. Because of the low emissivity of aluminum GCPs, they were visible in thermal images. Each GCP was measured with the stop-and-go technique through relative positioning by means of the NTRIP protocol (The Radio Technical Commission for Maritime Services, RTCM, for Networked Transfer via Internet Protocol) using two GNSS (global navigation satellite system) receivers. One of the receivers was a reference station for the GNSS Red Andaluza de Posicionamiento (RAP) network from the Institute for Statistics and Cartography of Andalusia, Spain, and the other, a Leica GS15 GNSS (Leica Geosystems AG, Heerbrugg, Switzerland), functioned as the rover receiver.  This weather station was equipped with a three-cup  anemometer, an air temperature and humidity sensor, and a barometer. Table 1 shows air temperature, percentage of relative humidity, mean wind speed, and atmospheric pressure during each UAV flight on this date. These UAV flights were used to acquire thermal images of soil and crop under different atmospheric and temperature conditions.

Thermal Image Processing
Thermal images were processed in three stages: thermal drift correction, geometric correction, and radiometric correction. Figure 3 describes the drift correction applied to the thermal images. As previously discussed, thermal drift is when the same location in the terrain presents different temperature values in different images. Thermal drift correction is based on points having a constant temperature during the UAV flight. In the first stage, a set of distinctive features are extracted from the UAV images using algorithms based on "structure from motion" (SfM) techniques described by Lowe [40]. SfM techniques extract individual features in each thermal image that are matched to their corresponding feature in the other images from the same UAV flight. As Figure 3 shows, a point appears in several images that belongs to different laps, each having a different temperature value due to drift effect. Every thermal image has a temporal reference, allowing the drift effect to be evaluated as it occurred in flight. In the proposed methodology, sensor drift is modeled as a function of time where each point of each image has a timestamp obtained by a GNSS sensor from the UAV autopilot. This methodology is applied to the set of characteristics extracted by the SfM algorithms. As a result, a mathematical model that describes thermal drift as it occurs for the duration of the UAV flight is achieved. Subsequently, this model is applied to all thermal images to obtain a new thermal image where temperature values are uniform on all corresponding points along the UAV flight. Six different drift correction models (DCMs) were developed to describe thermal drift for the duration of the UAV flight: exponential, exponential, lineal and polynomial order two, three, and four. Each DCM was applied to the UAV thermal images to generate a new collection of images, which were orthorectified and then processed into thermal orthomosaics. To obtain a single thermal orthomosaic of the area of interest, images have to be aerotriangulated, rectified, and finally mosaicked. Based on previous research results, Pix4dmapper by Pix4D SA was selected (Lausanne, Switzerland) and described by Mesas-Carrascosa [15] to do this processing. Afterward, digital numbers (DNs) of the thermal orthomosaics were converted to temperature values using the information obtained by the radiometric calibration of the sensor provided by the manufacturer. Finally, remotely sensed temperatures are influenced by the environmental conditions present at the time of UAV flight. Atmospheric correction of surface temperature is essential to extract absolute temperature measurements from thermal images, requiring the application of a radiative transfer model [41], a vicarious calibration [42], or an empirical method [43]. In this research, using an empirical method, two extreme temperature panels of 0.5 × 0.5 m were placed in the plot to record the hottest (black polymer panel) and the coldest (white polymer panel) temperatures on scene. Reference panels were close to five times larger than the GSD and, therefore, several homogeneous pixels appear in the thermal images. A temperature measurement of each reference panel was collected for each UAV flight with a Flir E60 heat gun (Flir Systems, Oregon, USA). These measurements were later used to correct the atmospheric effect on the thermal mosaics, applying an empirical line method, which defines a linear relationship between absolute temperature and sensor temperature [44].

Validation
The analysis of the proposed methodology was applied to all UAV flights both with and without thermal drift correction. A total of 37 georeferenced checkpoints (CPs) were placed along a transect perpendicular to the laps and crossing the middle zone of the plot and were read for temperature using a Flir E60 heat gun. These values were compared to the extracted values from the thermal orthomosaics obtained from the flights. The mean error and root mean square error (RMSE) were calculated for each model and flight. Moreover, Akaike's information criterion (AIC) [45] was used to identify the relative importance among all possible sets of DCMs per UAV flight where the best performing in each flight was identified by the lowest AIC score.
In addition, the correlation between flight time and thermal drift was defined. To do this, the thermal images in which a CP appeared were analyzed. Of all the possible thermal images where a CP appeared, the images with the most centrally situated CP were used in the mosaic process to obtain a thermal orthomosaic, as this is the standard method. As each image has an associated timestamp, each CP was temporally registered as it appeared in flight, which was then evaluated for the influence of flight time on thermal drift for each thermal orthomosaic obtained by a DCM. Figure 4 shows the variation of digital number per second (∆DN/t) along the duration of the UAV flights on features extracted from images and the different DCMs per flight developed for the proposal methodology. Comparing the four UAV flights, digital number per second in raw images varied both in absolute and relative terms. Therefore, the thermal sensor did not show a defined pattern in registering temperatures. Each DCM calculated per flight and its correlation coefficient (R 2 ) are summarized in Table 2. At 8:30 a.m. (Figure 4a), ∆DN/t had no defined behavior, showing a disperse distribution. Regardless of the mathematical model used, the correlation coefficient showed a value between 0.370 with the linear model and 0.382 for the exponential model of order 2. Therefore, there were no clear differences a priori between DCMs applied to thermal orthomosaics. At 12.30 p.m. (Figure 4b), ∆DN/t started to show a trend in its relation to flight duration as R 2 of some DCMs revealed. In this case, at the beginning of the UAV flight the ∆DN/t showed lower values over time (between 0 and 100 s) but then the variation increased, reaching the highest variation between 200 and 400 s. Thereafter, the variation decreased to values equal to the beginning of the UAV flight until 700 s where it again increased, although it did not reach the initial maximum values. In this flight, the bicubic model showed the highest R 2 , equal to 0.482, while the exponential model had the lowest, 0.190. At 16.00 p.m. (Figure 4c), ∆DN/t showed a defined evolution as the UAV flight progressed, which was reflected in higher R 2 values for all the DCMs. In this mission, ∆DN/ sec showed maximum values at the beginning of the UAV flight and then decreased over time. The exponential order 2 and lineal models had the lowest R 2 value (0.688) while the bicubic model had, again, the highest R 2 value (0.836). The other DCMs had an R 2 value higher than 0.7. Finally, the mission at 18:30 p.m. (Figure 4d) showed more oscillations during flight time; however, the range of ∆DN/ sec values was narrow. With this mission, maximum variations occurred at the beginning of the UAV flight, decreasing as the flight progressed. R 2 values were between 0.67 and 0.69 with no clear differences between the DCMs. Therefore, in analyzing the behavior of ∆DN/t along the UAV flight duration for these missions, it is possible to assert that the thermal sensor was not stable and that its operation varied in each UAV flight.

Results
Regarding R 2 , the DCMs obtained from the UAV flights at 16:00 p.m. and 18:30 p.m. had higher values than those from 8:30 a.m. and 12:30 p.m. because ∆DN/ sec showed a shorter range of values by time. This is due to the fact that ∆DN/t per second was different for each UAV flight, occurring with a lower frequency at later UAV flights.    The asterisks indicate the level of significance (* p < 0.05, ** p < 0.001, n.s. not significant). Figure 5 shows the thermal orthomosaic histograms for each flight mission with the applied DCMs, as well as without a DCM; Table 3 summarizes the statistics for each. First, the histograms manifest a clear difference between the thermal orthomosaics where a DCM was applied compared to those where no DCM was applied. Moreover, the temperature distribution was quite similar in all the thermal orthomosaics where a DCM was applied. Because the study area had two differentiated classes, vegetation and bare soil, a bimodal distribution was expected to describe temperature distribution of the scene. However, at the 8:30 a.m. UAV mission (Figure 5a), all thermal orthomosaics had a normal distribution as Sarle's bimodality coefficient (SBC) showed. The temperature range on those thermal orthomosaics where a DCM was applied ranged from 15 to 35 • C, 20 • C occurring with the most frequency. Conversely, the thermal orthomosaic without any drift correction showed a broader temperature range from 15 to 43 • C and a right-skewed distribution. Instead of a bimodal distribution, a normal distribution of temperatures occurred in the early morning on the thermal orthomosaics with drift correction because the bare soil and vegetation had not yet absorbed heat from the sunlight and consequently the temperatures of both were similar. Therefore, at this hour, both classes showed no clear difference in thermal behavior. Conversely, at 12:30 p.m. (Figure 5b), the UAV flight histograms showed two different shapes irrespective of whether a DCM was applied or not. However, when no drift correction was applied, the temperature distribution was, again, similar to a normal distribution with SBC being less than 5/9. As a result, the non-corrected histogram did not properly mark bare soil and vegetation with different temperatures. On the other hand, all drift-corrected thermal orthomosaics showed an SBC higher than 5/9, having a bimodal distribution with two differentiated peaks. In this mission, bare soil and vegetation had different behaviors as they correlated to sunlight. Although both classes increased their temperature, vegetation (left peak) showed a mean temperature equal to 30 • C, which was lower than bare soil (right peak), which reached a mean value equal to 50 • C. At 16:30 p.m. (Figure 5c), both vegetation and bare soil increased in temperature, which was properly shown when a DCM was applied. If corrections were not applied, the temperatures were also higher but the classes were not accurately represented in the histogram.
A comparison of the histograms shows that vegetation temperature increased as the day progressed, reaching the highest temperature at 18:30 p.m. while soil temperature increased until 16:30 p.m. and then began to decrease. At 18:30 p.m. (Figure 5d), the left peak of the histogram has a higher frequency than the right peak as the vegetation maintained a stable temperature while the bare soil temperature decreased, which also caused the distance between the peaks to be reduced. These results are because the soil was cooling due to the declining sun and greater shadow cover.
All of the DCMs used successfully described this occurrence from 12:30 p.m. to 18:30 p.m. Moreover, the histograms of each DCM for every flight had a similar distribution with the quartic model presenting the greatest differences in portions of the flights. The histograms show that from 12:30 p.m., the temperature difference between the vegetation and bare soil was quite distinctive and remained so until 16:30 p.m. when the difference began to decrease. Therefore, as in Bellvert et al. (2014) [43], it is recommended that thermal UAV flights with agronomic objectives are performed between 12:00 and 16:30 p.m. However, without the method proposed in this paper, non-corrected thermal orthomosaics will not sufficiently differentiate soil and vegetation.

Validation
Once the histograms of the drift corrected thermal orthomosaics accurately described the presence of vegetation and bare soil in the study area, the next step was to analyze which DCM had the greatest thermal accuracy. Figure 6 illustrates the 16:30 p.m. thermal orthomosaics with the applied DCMs, which are detailed in Table 2, as well as without drift correction. When drift correction was not applied (Figure 6g), the resulting thermal orthomosaics ordered the temperature values. From a visual analysis, neither bare soil nor vegetation showed stable temperature. Moreover, temperature changed along the north and south laps, registering higher temperature values as the UAV flight progressed. This effect is pronounced in this paper because NUC was switched off to obtain an extreme example; in other cases, the drift effect would be less pronounced. This also explains the pronounced skewness in the non-corrected histograms presented in Figure 5. Based on the authors' results, the temperature variation for this sensor under normal conditions was less than 0.5 • C per minute, similar to the variations reported by Olbrycht and Więcek (2015) [29].
Regarding the thermal orthomosaics where a DCM was applied (Figure 6a-f), no visual temperature differences from north to south were detected. Instead, the temperature variations were linked to the state of vegetation and bare soil as explained by the histogram analysis above. In addition, comparing all of the thermal orthomosaics generated with a DCM resulted in similar orthomosaics with the exception of the quartic model (Figure 6f). This model generated colder temperatures in the south of the plot compared to the north of the plot, which was not present in the other DCMs. These differences were not detected in the field campaign, suggesting that the quartic model did not adequately describe thermal drift in the UAV flights. This result occurred for all of the UAV flights assessed. Table 4 summarizes the results of thermal quality control on the thermal orthomosaics from each UAV flight with and without the applied DCMs by mean error and standard deviation (SD) and Akaike's information criterion (AIC). In addition, a correlation coefficient (r 2 ) between error and flight progress was calculated. In all the orthomosaics, as expected, where a DCM was not applied, the temperature errors were higher than where a DCM was applied. Therefore, it is necessary to pre-process thermal images taking into account the behavior of the microbolometer registering temperature values. The 8:30 a.m. mission showed a higher rate of error than the other UAV missions independent of which DCM was applied. At this time, the error ranged from 0.88 ± 0.8 • C using the bicubic model correction to 1.01 ± 0.81 • C using the exponential correction model. These higher errors are explained by the different environmental conditions while performing the UAV flight and measuring ground truth. In the early morning, the sun is ascending and an object's or coverage's superficial temperature changes in a short time interval. Although ground truth measurements were performed immediately after the UAV flight finished (twenty minutes), the duration of the UAV flight was long enough for temperature values of the olive trees and bare soil to change both for images and for on-ground measurements, influencing this mission's results. At 12:30 p.m., the mean error decreased to values around 0.2 • C ± 0.5 • C for all DCMs applied with the cubic model generating the smallest error (0.10 • C ± 0.45 • C) while the exponential model order 2 had the highest error (0.29 • C ± 0.56 • C). Conversely, at 16:30 p.m., the errors differed depending on the DCM used. For this mission, the quartic model showed the highest error (1.57 • ± 1.22 • C) while the cubic model had the lowest (0.06 ± 0.45 • C). Finally, at 18:30 p.m., the cubic model again showed the greatest accuracy with an error equal to 0.26 ± 0.58 • C and the quartic model being the worst with an error equal to 1.55 ± 0.85 • C. Even so, the highest mean error and SD were found in those thermal orthomosaics where no DCM was applied, independent of the mission. Moreover, the calculated AIC scores identified the cubic model as the most consistent DCM, showing the minimum score in each UAV flight.  To demonstrate, Figure 7 shows the relationship between temperature obtained from the UAV flight and the on-ground measurements at 16:30 p.m. The continuous line represents a perfect correlation of temperatures between both sets of temperatures and a discontinuous line represents the adjusted lineal model from both sets. When no DCM was applied (Figure 7g), the temperature values in the thermal orthomosaic did not show any relationship with those on the ground having a broad range of variation. On the other hand, in general, when a DCM was applied, temperature values on the thermal orthomosaics and their corresponding on-ground measurement were similar. However, there were differences between the DCMs. The quartic model (Figure 7f) yielded the highest deviation from a perfect correlation although not as broad as in the case of not using a DCM. The same occurred when the exponential model (Figure 7a) or lineal model (Figure 7c) was applied, although not as evident as in the previous case. The remaining DCMs considered did not show significant differences in their temperature values. Moreover, to have an acceptable range of temperature error, it is necessary that this error occurs independently of the flight duration, meaning that the drift effect of the microbolometer has been adjusted accordingly. Table 4 shows this relationship through a correlation coefficient and its Pearson analysis and Figure 8 shows the relation between error and flight duration for each thermal orthomosaic obtained from the 16:30 p.m. mission where a DCM was applied. As expected, the thermal orthomosaics where a DCM was not applied showed high R 2 between error and flight duration in all UAV flights, meaning that the error is dependent on flight duration. As to the thermal orthomosaics where a DCM was applied, the quartic model was the one whose errors showed a significant correlation with flight duration in all UAV flights. The other DCMs had no significant correlation at 8:30 a.m., 12:30 p.m., and 18:30 p.m. At 16:30 p.m., the cubic model was the only DCM whose results were independent of flight duration. Therefore, when considering the relationship between temperature error and flight duration, the cubic model adequately described drift effect on temperature measurements on all performed UAV flights. In this study and our experience, using the Gobi 640 thermal sensor on subsequent UAV flights has shown that the cubic drift model is the most adequate. However, for other thermal sensors, it is recommendable to analyze which model would better describe drift effect using the proposed methodology. The temperature error obtained in this research when a DCM was applied is equal to those described in the literature [8,38] and therefore presents itself as another option for agronomical projects. Moreover, the proposed methodology presents an advantage as it optimizes flight time. In other studies, it is necessary to stabilize the sensor in flight [36] or it is necessary to recurrently fly over temperature targets for reference to calibrate thermal images [38]. These two strategies spend the limited battery charge and, therefore, reduce the area covered by the UAV flight due to the decrease in useful flight time. With the proposed methodology, flight time is maximized without loss of accuracy in temperature values. Moreover, it allows the use of a thermal sensor on board the UAV regardless of the drift effect.
The UAV flights were performed under stable weather conditions and were of short duration (20 min) and as such, there were no changes in the surface temperature. However, UAV flights with longer duration (longer than 30 min) and/or under unstable weather conditions must be evaluated in future work due to changes in atmospheric or sunlight conditions. In these cases, in addition to the drift correction, it would be necessary to have an instantaneous atmospheric correction adapted to the UAV flight conditions. One possible methodology is to equip the UAV platform with additional sensors that record values of temperature, relative humidity, atmospheric pressure, incident radiation, and wind speed. These parameters would allow the possibility of determining the atmospheric conditions linked to each individual thermal image and, therefore, at each moment of flight time. These parameters along with the drift effect model would allow more precise and accurate temperature values.

Conclusions
Remote sensing using lightweight uncooled thermal sensors on board UAVs is a useful tool for measuring crop temperatures. However, drift effect on registered temperatures can invalidate its agronomical applications where an accuracy greater than 1 • C is necessary. The present research has developed methodology to remove drift effect on temperature using a lightweight microbolometer thermal sensor on board a UAV. In this study, removing drift effect on thermal images is based on redundant information around objects that appear in different overlapping images from a UAV flight that covers the area of interest. Different mathematical models were explored to describe drift effect with the cubic drift model yielding the best results on separate missions performed for this research.
These models were tested in four UAV missions at different hours at the same location. If no drift correction was applied, the thermal orthomosaics did not adequately describe crop temperatures, invalidating their use within an agronomical context. Contrarily, if a drift correction model was applied using the proposed methodology, the results improved with a range of error that would be adequate for agronomical projects. Moreover, the accuracy is in the same range as other authors' results but with the added benefit of not requiring any special in air UAV flight operation and thusly increasing useful flight time. This is an important point to those UASs with short flight durations due to limited battery power. In addition, this methodology was applied to a single UAV flight and, therefore, the proposed drift correction model is adaptable to specific flight conditions.
In this study, the cubic drift model offered the best results. However, the authors recommend exploring the behavior of a particular uncooled thermal sensor to determine which model would best describe the drift effect using the proposed methodology.