Self-calibration of UAV Thermal Imagery Using Gradient Descent Algorithm

1. Introduction

1.1. Uncooled thermal cameras issues

Thermal imagery from UAVs has several advantages in environmental applications as they can provide a low-cost, high-resolution, and flexible approach to environmental monitoring and management. UAV-based thermal and narrowband multispectral imaging sensors can provide low-cost approaches to meet the critical requirements of spatial, spectral, and temporal resolutions for vegetation monitoring [1]. UAVs provide high spatial resolution and flexibility in acquisition and sensor integration, which can be used for land cover classification, change detection, and thematic mapping [2]. Aerial thermography from low-cost UAVs can be used to generate digital thermographic digital terrain models, which find application in the classification of land uses according to their thermal response [3]. Thermal remote sensing has many potential uses in precision agriculture, including monitoring plant hydration levels, identifying instances of plant diseases, evaluating crop yield, and analysing plant characteristics [4].

Limitations related to the maximum load of the UAV, as well as cost constraints, lead to the use of mainly lightweight uncooled thermal cameras. This solution has several shortcomings when measuring temperature under field conditions. Uncooled UAV thermal cameras require careful calibration and correction for various factors to obtain accurate temperature measurements. They suffer from vignette effects, sensor drift, ambient temperature influences, and measurement bias, which can be corrected with an ambient temperature-dependent radiometric calibration function [5]. Non-radiometric uncooled thermal cameras are highly sensitive to changes in their internal temperature and require empirical line calibration to convert camera digital numbers to temperature values [6]. Fluctuations in the temperature of the focal plane array (FPA) array detector, wind, and irradiance can affect temperature measurements, and that adequate settings of camera gain and offset are crucial to obtaining reliable results [7]. Uncooled thermal cameras also suffer from thermal-drift-induced nonuniformity or vignetting [8]. The above factors contribute to fluctuations in camera’s accuracy from ±0.5 °C in laboratory conditions to ±5 °C in unstable UAV flight conditions [6].

The literature review provides a number of solutions to reduce this problem. To minimise measurement errors in UAV thermal cameras, it is suggested that the camera be warmed up for 15-40 minutes before starting the actual measurement [6,8]. There are also several methods to eliminate unwanted phenomena outcomes with calibration under controlled conditions. Ribeiro-Gomes et al. proposed a calibration algorithm based on neural networks that allows for increased measurement accuracy [9]. In another study by Yuan and Hua, nonuniformity of vignetting was examined in UAV-based uncooled thermal cameras and a simple vignetting reduction method using reference images of a homogeneous target was proposed [8]. In Aragon et al., an ambient temperature-dependent radiometric calibration function was used to correct for sensor drift, ambient temperature influences, measurement bias, and vignette effects [5]. An advanced radiometric calibration approach that stabilises the camera's response, removes fixed pattern noise, and converts thermal image values into object temperatures was proposed by Lin et al. [10]. The disadvantage of calibration-based methods is that they require non-standard equipment and the extra effort needed to collect calibration data. A correction based only on the field data collected by the UAV system (self-calibration) would be favoured. Mesas-Carrascosa et al. developed a method for bias correction using redundant data from overlapping areas of aerial thermal images using mathematical modelling of the metric of “variation of digital number per second” [11]. Unfortunately, despite good results, the authors did not provide a specific explanation, a formula, or sample source code illustrating the calculation of this metric. Also, since this method leverages the rate of digital number variation, it requires precise timing of image acquisition achieved with a custom drone attachment. Often by default the time information available in the image metadata is provided with a precision of a second, which is not sufficient for rate of digital number variation calculation, as during the flight successive photos are taken even as often as every ca. 1 second.

Another important aspect of the problem is the choice of thermal imaging camera brand. Some manufacturers of cameras prevail in publications related to the topic discussed. The files generated by such cameras offer a great deal of freedom to recalibrate temperature measurements, as the sensor response raw data can be easily accessed. Another manufacturer offering in recent years its own thermal imaging cameras dedicated to drones, has another internal processing algorithms policy. In their case, raw sensor data is not openly available, and it is impossible for the user to recalibrate the device. The user has to rely on producer's closed-source Thermal SDK library allowing to make corrections to the measurement taking into account air humidity, target emissivity or flight altitude. It is not known what exact processes are represented in this algorithm. Moreover, we noticed its significant shortcomings: i) the maximum flight altitude acceptable by the library is 25m, which is too low for most flights covering large areas, ii) the library returns errors for some humidity values (in our tests, most often for relative humidity of about 60%), iii) for a given frame it is possible to select only one emissivity (it is not possible to set different emissivities for different objects in the same photo). However, most of the UAV dedicated cameras are susceptible to adverse phenomena typical of uncooled thermal cameras (bias, vignette effect).

1.2. Gradient descent algorithm

A major source of measurement uncertainty is bias, which can be addressed by investigating redundant data from overlapping areas of the images. It is desirable that the temperature difference between all pairs of overlapping images be as small as possible. This is a minimisation problem that can be solved with any optimisation method. This work focusses on the use of gradient descent, which is an optimization algorithm commonly used in machine learning. However, beyond its main application in machine learning, gradient descent can be used to optimise any differentiable objective function.

Gradient descent works by iteratively updating a set of parameters in the direction of the steepest descent of a cost function. The algorithm computes the gradient of the cost function with respect to the parameters, and then updates the parameters by taking a step in the direction of the negative gradient. The learning rate parameter determines the size of the step and it is usually set to a small value to avoid overshooting the minimum of the cost function. The algorithm continues to iterate until the cost function converges to a minimum or a stopping criterion is met. With the popularity of machine learning, gradient descent has become more accessible thanks to the development of software libraries that accelerate the algorithm by using the GPU (Graphical Processing Unit).

1.3. Objective of this work

The aim of this work is to develop a method for post-processing of aerial thermal imagery that reduces the effects of undesirable phenomena occurring in uncooled thermal cameras without the use of non-standard equipment, such as reference black bodies or custom UAV attachments, and without access to raw thermal sensor data. Achieving this goal is important, as the use of thermal imaging cameras is becoming widespread, but proven post-processing methods are lacking.

2. Materials and Methods

2.1. Data

The thermal pictures were collected by the DJI Matrice 300 RTK drone system equipped with a Zenmuse H20T multicamera sensor that contains the thermal camera. Data used in this study was collected of several areas in southern Poland:

• Area around ca. 500 m Kocinka stream stretch near Grodzisko village (50.8715 N, 18.9661 E)
• Area around ca. 350 m Kocinka stream stretch near Rybna village (50.9371 N, 19.1134 E)
• Area around ca. 160 m Sudół stream stretch near Kraków city (50.0999 N, 19.9027 E)

The water temperature of the rivers was measured directly along their course in each case. The results were constant for the entire section and did not depend on chainage. Table 1 provides details of the locations, dates, conditions of the surveys, and measured temperatures.

2.2. Algorithm

2.2.1. Vignette correction

The correction of the vignette effect was conducted using the “single image” method [12]. ‘Single image’ means that the algorithm tries to model the correction of the vignette effect only on the basis of the currently processed image and does not have any auxiliary data available (e.g., an image of a homogeneous target with a clearly visible vignette effect). It was implemented by translating the MATLAB code available at https://github.com/GUOYI1/Vignetting_corrector into Python. Several changes have been made to adapt the algorithm to work with thermal images. Standard images are represented with floating point values from 0 to 1, or integers from 0 to 255. The original implementation of the algorithm assumed this data format, so it had to be modified to work on unconstrained float values. Thermal image values can also be negative. The vignette correction algorithm works on logarithms of pixel values, so temperatures are converted to Kelvins in order to avoid the logarithm of 0. The vignette correction algorithm tends to increase the brightness of the photo. For standard photos this is not a problem, but for thermal images, where increasing the brightness will result in a bias in the temperature reading, this cannot be accepted. We assume that images with vignette effect occurring contain correct temperature in their central part. Therefore, we can use the central part of the photo before correction as a reference level for the photo obtained after applying the vignette correction algorithm. The bias is compensated for by subtracting the average difference between the central areas of the images before and after vignette correction. While the vignette effect mainly affects the edge part of the image, the central area size is defined as a rectangle with sides twice smaller than the whole image.

2.2.2. Georeferencing

The georeferencing algorithm is largely based on the aerial photos stitching solution proposed by Luna Yue Huang [13]. Modifications in our method consists of direct georeferencing to Universal Transverse Mercator (UTM) coordinate system.

The first part of the georeferencing procedure is initial georeferencing of images based on EXIF metadata. These contain information about the geographic coordinates, yaw, and altitude of the drone camera at the time the image was taken. This information, along with the angle of field of view taken from the camera specifications, allows estimation of an affine transformation parameters of translation ($v_x$, $v_y$), scale ($s_x$, $s_y$), and counter clockwise rotation ($\theta$) that allows to embed the image in an UTM coordinate system. This is not an accurate estimation, due to the low accuracy of the UAV's GPS geolocation and high sensitivity to the shift of the camera viewing area under the influence of wind blows. In this estimation the scale $s_x$ and $s_y$ parameters are equal, and the shear ($c_x$, $c_y$) coefficients are 0. The transformation parameters allow us to obtain the transformation matrix $A$ based on Equation 1.

$$\begin{gathered} A=\left(\begin{array}{ccc}\cos \theta & -\sin \theta & v_x\\ \sin \theta & \cos \theta & v_y\\ 0& 0& 1\end{array}\right)·\left(\begin{array}{ccc}{s}_x& c_x& 0\\ c_y& s_y& 0\\ 0& 0& 1\end{array}\right) \\ =\left(\begin{array}{ccc}{s}_x\cos \theta -c_y\sin \theta & c_x\cos \theta -s_y\sin \theta & v_x\\ c_y\cos \theta +s_x\sin \theta & s_y\cos \theta -c_x\sin \theta & v_y\\ 0& 0& 1\end{array}\right), \end{gathered}$$

(1)

Based on the initial georeferencing, pairs of overlapping images are found. To ensure that all possible pairs are found even when errors in the initial georeferencing may indicate a lack of overlapping, footprints are expanded (buffered) by a dilation operation by the given padding value. Each pair of images is aligned and a relative transformation matrix $A_R$ is found by a well-established in the field of image vision stitching method: i) finding corresponding points (keypoints) in both images using scale-invariant feature transform (SIFT) [14], ii) matching keypoints using fast approximate nearest neighbour search (FLANN) [15], and iii) estimation of best transformation matrix using random sample consensus (RANSAC) [16]. Relative transformations are verified in two stages: i) since we assume that the images are taken from the same altitude, the scaling factor in the relative transformation of two overlapping images must be approximately equal to 1. If the value of the scaling factor for relative transformations is outside the range 0.9 – 1.1, the relative transformation is considered incorrect. ii) as all images are taken in nadir orientation, the relative transformations should not perform a shearing operation. When there is no shear, the absolute values of the relative transformation matrix $A_R$ satisfy the following condition: $\left|A_{R11}\right|=\left|A_{R22}\right|$ and $\left|A_{R12}\right|=\left|A_{R21}\right|$. If deviation between these values is greater than 0.1, the relative transformation is considered to be incorrect. If the relative transformation is found to be incorrect, according to Luna Yue Huang solution, an attempt is made to establish a relative transformation using the same pair of images, but resized to twice lower resolution. Such a procedure allows you to obtain a different keypoints using for scale determination. If the relative transformation obtained in this way again fails to pass verification, the pair of photos is discarded.

Based on the obtained pairs, an undirected graph was constructed [17], where each node is an aerial photo, and each edge is a valid relative transformation between a pair of photos. The connected components are then extracted from the graph. A connected component in graph theory is a group of vertices in a graph that are all directly or indirectly connected to each other. Only images (nodes) from a connected component containing the largest number of nodes are used for further processing.

Iterating over each pair of images, the coordinates of the corners of the second image are calculated in the pixel coordinate system of the first image using equation 2. This idea is also visualised in Figure 1.

$$p_{1st}=A_R·p_{2nd},$$

(2)

where:

$A_R$ – relative transformation matrix between 2nd image and 1st image pixel coordinate system,

$p_{2nd}$ – coordinates of corner point expressed in 2nd image pixel coordinate system,

$p_{1st}$ – coordinates of the same point expressed in 1st image pixel coordinate system.

Optimisation of georeferencing of all images involves tuning the absolute geographic transformation parameters ($v_x$, $v_y$, $s_x$, $s_y$, $\theta$) to make the relative transformations recovered from them as close as possible to the relative transformations obtained earlier using alignment of each pair separately. Recovery of the relative transformation matrix $\widehat{A_R}$ from the absolute geographic transformation matrices $\widehat{A_{1st}}$ and $\widehat{A_{2nd}}$ of two georeferenced images is obtained according to Equation 3.

$$\widehat{A_R}={\widehat{A_{1st}}}^{-1}·\widehat{A_{2nd}},$$

(3)

$c_x$ and $c_y$ shear parameters are not tuned during the optimization. Although tuning of these parameters further improves the matching between pairs of photos, it also introduces an unacceptable shearing of the entire mosaic of photos.

Similarly to equation 2, the coordinates of points expressed in the pixel reference system of the 2nd photo can be converted to coordinates expressed in the reference system of the 1st image using equation 4.

$$\widehat{p_{1st}}=\widehat{A_R}·p_{2nd},$$

(4)

where:

$\widehat{A_R}$ – relative transformation matrix between 2nd image and 1st image pixel coordinate system,

$p_{2nd}$ – coordinates of point expressed in 2nd image pixel coordinate system,

$\widehat{p_{1st}}$ – coordinates of the same point expressed in 1st image pixel coordinate system.

Optimisation is performed using the gradient descent method with the loss function $\mathcal{L}$ that consist of two components ${\mathcal{L}}_R$ and ${\mathcal{L}}_A$. Component ${\mathcal{L}}_R$ is the mean Euclidean distance between the points at the corners of the images obtained from the relative transformations recovered from the tuned absolute geographic transformations and the points from the corners of the images from the relative transformations obtained during the pair alignment.

$${\mathcal{L}}_R=\frac{1}{N_P}\sum _{i=1}^{N_P}\frac{1}{4}\sum _{j=1}^4‖p_{1st,i,j}-\widehat{p_{1st, i,j}}‖,$$

(5)

where:

$N_P$ – number of pairs,

$p_{1st,i,j}$ – point of j-th corner of the i-th pair 2nd image expressed in the 1st image pixel coordinate system estimated from pairs alignment,

$\widehat{p_{1st,i,j}}$ – point of j-th corner of the i-th pair 2nd image expressed in the 1st image pixel coordinate system estimated using relative transformation recovered from absolute geographic transformations tuned during optimization.

Component ${\mathcal{L}}_A$ is a mean Euclidean distance between geographic centroids of tuned image footprints and geographic centroids of image footprints obtained during initial georeferencing using EXIF data.

$${\mathcal{L}}_A=\frac{1}{N_I}\sum _{i=1}^{N_I}‖p_i-\widehat{p_i}‖,$$

(6)

where:

$N_I$ – number of images,

$p_i$ – point of centroid of i-th image obtained from EXIF data expressed in geographic coordinate system,

$\widehat{p_i}$ – point of centroid of i-th image obtained from absolute geographic transformation tuned during optimization.

By minimizing the ${\mathcal{L}}_R$ factor, the optimized absolute geographic transformations are adjusted so that the retrieved from them relative transformations between pairs of images are as close as possible to the relative transformations previously obtained separately for each pair in the alignment process. Minimizing the ${\mathcal{L}}_a$ factor ensures that the entire mosaic will not move in any direction during the optimization. Equation 7 is the final formula for the $\mathcal{L}$ loss function. To minimise the impact on the relative matching of images using ${\mathcal{L}}_R$ component, the ${\mathcal{L}}_A$ component is multiplied by a factor of ${10}^{-6}$. The factor value was selected experimentally by qualitative assessment of results tested for different values.

$$\mathcal{L}={\mathcal{L}}_R+{10}^{-6}{\mathcal{L}}_A,$$

(7)

2.2.3. Thermal calibration

In order to minimise the temperature bias between overlapping images using Gradient Descent, a consistent dataset has to be prepared. If the reduction of the vignette effect has not yielded sufficient results and the photos overlap is large enough, the areas of the photos near the edges where the vignette effect affects the temperature images the most can be truncated (Figure 2).

For the two images from each pair, only the common part is retained, as well as a mask that allows us to reproduce the irregular shape of the cutout from the rectangular array. Temperature cutouts of both images and mask are resized to array of $32\times 32$ pixels. During resizing, the temperatures are interpolated using the bilinear method, and the mask is interpolated using the nearest neighbours method. The example result of the common clip is shown at Figure 3. The example corresponds to the images shown in Figure 1 and Figure 2. The rotation is due to the transformation to an absolute geographic datum. The lack of the preservation of aspect ratio is due to the need for scaling to fill the whole area of square array of $32\times 32$.

By resizing the slices to $32\times 32$ pixels, it was possible to create a coherent dataset of three arrays (masks and pair of images common part) of size $N_p\times 32\times 32$ ($N_p$ is number of pairs), which could then be used in optimization using the gradient method.

Optimisation of the temperature values involves tuning the $b$ parameter that is used for the offset of each image according to equation 8.

$$\widehat{v}=v+b,$$

(8)

where:

$\widehat{v}$ – temperature value with applied calibration,

$v$ – temperature value before calibration,

Optimisation uses a loss function $\mathcal{L}$ that consists of two components: ${\mathcal{L}}_B$ and ${\mathcal{L}}_M$. Component $L_B$ is pixelwise mean squared error between calibrated temperature values of first and second cutout images from given pair. It is calculated with Equation 8.

$${\mathcal{L}}_B=\frac{1}{N_P}\sum _{i=1}^{N_P}\frac{1}{32}\sum _{j=1}^{32}\frac{1}{32}\sum _{k=1}^{32}{\left(\widehat{v_{1st, i,j,k}}-\widehat{v_{2nd, i,j,k}}\right)}^2,$$

(9)

where:

$N_P$ – number of pairs,

$\widehat{v_{1st, i,j,k}}$ – value of calibrated first image in i-th pair at pixel coordinate $(j,k)$,

$\widehat{v_{2nd, i,j,k}}$ – value of calibrated first image in i-th pair at pixel coordinate $(j,k)$,

${\mathcal{L}}_M$ is absolute difference between mean of all uncalibrated images and mean of calibrated images. It is calculated with Equation 9.

$${\mathcal{L}}_M=\left|\mu -\widehat{\mu }\right|,$$

(10)

where:

$\mu$ – mean value of uncalibrated images,

$\widehat{\mu }$ – mean value of calibrated images.

By minimizing the ${\mathcal{L}}_B$ component, the bias between the temperature values in the overlapping images is minimized. The ${\mathcal{L}}_M$ component ensure that during optimization the whole mosaic will not drift from initial mean value. The total loss is calculated using equation 11. To minimise the impact on bias correction using the ${\mathcal{L}}_B$ component, the ${\mathcal{L}}_M$ component was multiplied by factor of ${10}^{-2}$. The factor value was selected experimentally by qualitative assessment of results tested for different values.

$$\mathcal{L}={\mathcal{L}}_B+{10}^{-2}{\mathcal{L}}_M,$$

(11)

2.3. Final processing and result analysis

The rivers in the measured areas have homogenous temperatures throughout the surveyed sections and a well-defined emissivity value of 0.95. Therefore, they can be used to verify the proposed method. The sets of both the original and bias-corrected images were offset (the same offset value for all images from the given set) so that the average water body temperature read from the thermal mosaic was equal to the measured reference water temperature. We decided on this because the temperatures returned by the DJI Thermal SDK may differ significantly from the actual temperatures, it is straightforward correction and makes it easier to compare measurements from before and after calibration.

3. Results

3.1. Vignette correction

A 'single image' algorithm was used to correct the vignette effect. For this reason, it is prone to misinterpretations of the image and the resulting erroneous correction of the vignette effect. Figure 4 shows an example of correct vignette effect correction. In this case, the standard deviation of the image as a result of the vignette effect correction has decreased, which manifests itself as a narrowing of the histogram of pixel values. Figure 5 shows an incorrect correction of the vignette effect, where it is likely that the river near the left edge of the image caused the misinterpretation. In this case, the correction resulted in an increase in the standard deviation of the photo, which manifests itself as a widening of the histogram of pixel values.

To quantitatively investigate the effectiveness of the devignetting algorithm, it was checked for how many images a reduction in standard deviation was obtained as a result of devignetting. For 74.4% of the images (2261 out of 3037 total images) there was a reduction in the standard deviation. For 25.6% of the images (776 out of 3037 total images) there was an increase in the standard deviation. The average reduction in standard deviation was -0.07 °C while the average increase in standard deviation was 0.01 °C. The detailed distribution of the changes in standard deviation as a result of the vignette effect correction is shown in the histogram in Figure 6.

3.2. Visual assessment

Figure 7, Figure 8, Figure 9, Figure 10 and Figure 11 show parts of thermal image mosaics before and after calibration for each case study. These are not orthomosaics - they were created straightforwardly by displaying all georeferenced images. Where the images overlap, a later image from the collection is displayed.

3.3. Waterbody temperature

Table 2 summarizes the obtained RMSE and MAE of stream temperature measurements using uncalibrated and calibrated thermal images. Both calibrated and uncalibrated images were subjected to the offset described in section 2.3.

Figure 12 show the temperatures sampled along the river centreline from the uncalibrated and the calibrated images. There are multiple measurements for a single chainage value, as a given point may be present in multiple overlapping images. Peaks in the temperature readings that have not been compensated by calibration are sampled on vegetation overhanging the water surface.

4. Discussion

4.1. Devignetting

The single-image vignette correction method has proven successful when applied to thermal images. Although there was an increase in the standard deviation for about a quarter of the images. This was a minor increase, as on average a seven times greater reduction in standard deviation was recorded for the rest of the three-quarters of the images.

4.2. Visual assessment of mosaics

A visual comparison of the mosaics before and after calibration shows an improvement in temperature fluctuations, particularly between successive flight passages. The algorithm also handles the correction of undervalued images taken at the beginning of the flight, when the camera is warming up. Correction of the images from the case study 20211215_kocinka_rybna and 20220118_kocinka_rybna has allowed us to see the warmer tributary of the stream, which is fed by warmer groundwater (Figure 7 and Figure 8 after calibration).

4.3. Measurement precision and accuracy

The proposed solution based on minimising the bias between overlapping images significantly increased the precision of the water temperature measurements taken in all case studies. Despite the improvement in calibration precision, the accuracy can still be poor due to imperfections of the Thermal SDK library provided by the manufacturer, which was used for conversion of raw data to temperature values. However, an increase in accuracy can be easily achieved by an additional offset applied to all images (see Section 2.3) if any reference points are available for the study area. After applying this offset, the average root mean squared error was reduced by 39.0% and the mean absolute error by 40.5%.

4.4. Applications

Previous methods of reducing adverse phenomena occurring in uncooled thermal cameras have relied on calibration under controlled laboratory conditions and required non-standard instrumentation such as reference black bodies or custom UAV attachments. For this reason, they require an additional effort of specialised operator and cannot be automated. On the other hand, our proposed solution requires minimal effort as it only requires data collected during a standard UAV flight and is highly automated making it possible to use by a non-specialised operator. This also makes the solution ready to be implemented as part of a photogrammetric software workflow.