Disaggregation of Landsat-8 Thermal Data Using Guided SWIR Imagery on the Scene of a Wildfire

Abstract

Thermal data products derived from remotely sensed data play significant roles as key parameters for biophysical phenomena. However, a trade-off between spatial and spectral resolutions has existed in thermal infrared (TIR) remote sensing systems, with the end product being the limited resolution of the TIR sensor. In order to treat this problem, various disaggregation methods of TIR data, based on the indices from visible and near-infrared (VNIR), have been developed to sharpen the coarser spatial resolution of TIR data. Although these methods were reported to exhibit sufficient performance in each study, preservation of thermal variation in the original TIR data is still difficult, especially in fire areas due to the distortion of the VNIR reflectance by the impact of smoke. To solve this issue, this study proposes an efficient and improved disaggregation algorithm of TIR imagery on wildfire areas using guided shortwave infrared (SWIR) band imagery via a guided image filter (GF). Radiometric characteristics of SWIR wavelengths could preserve spatially high frequency temperature components in flaming combustion, and the GF preserved thermal variation of the original TIR data in the disaggregated result. The proposed algorithm was evaluated using Landsat-8 operational land imager (OLI) and thermal infrared sensor (TIRS) images on wildfire areas, and compared with other algorithms based on a vegetation index (VI) originating from VNIR. In quantitative analysis, the proposed disaggregation method yielded the best values of root mean square error (RMSE), mean absolute error (MAE), correlation coefficient (CC), erreur relative globale adimensionelle de synthèse (ERGAS), and universal image quality index (UIQI). Furthermore, unlike in other methods, the disaggregated temperature map in the proposed method reflected the thermal variation of wildfire in visual analysis. The experimental results showed that the proposed algorithm was successfully applied to the TIR data, especially to wildfire areas in terms of quantitative and visual assessments.

1. Introduction

Thermal data is related to emissivity of objects on the ground, and this property can account for biophysical phenomena [1,2]. Recently, thermal remote sensing based on satellite sensors has made it possible to establish temperature maps over wide areas. However, a trade-off exists between spatial and spectral resolutions in any spaceborne remote sensing system [3,4]. This limitation appears especially in thermal infrared (TIR) sensors. For instance, TIR sensors in the Landsat-8 platform obtain the radiation emitted by the earth in a range of 10–13 μm, but the TIR band retrieves lower total radiation than visible-near infrared (VNIR) and shortwave infrared (SWIR) bands. Due to this radiometric characteristic of the TIR sensor, the spatial resolution of the TIR band is coarser than that of other bands. In addition, the limited spatial resolution of TIR sensors emphasizes the necessity for disaggregation of TIR imagery, also called thermal sharpening or downscaling, in order to improve the spatial resolution of the thermal data products [5]. The purpose of TIR disaggregation is to sharpen the coarser spatial resolution of TIR images by introducing spatial detail from a finer resolution multispectral (MS) image, without distortion of spectral information.

To achieve this purpose, various disaggregation methods of TIR imagery have been proposed in the last few decades. Generally, most have been focused on the relationship between low-resolution TIR and high-resolution vegetation index (VI), especially normalized difference vegetation index (NDVI), because densely vegetated areas tend to have relatively low surface temperature [6,7,8,9,10,11,12,13,14,15]. Disaggregation procedures for radiometric surface temperature (DisTrad) [6] and thermal sharpening (TsHARP) [7] algorithms are viewed as the leading methods due to their effectiveness and simplicity. Based on the TsHARP algorithm, one category of the refined disaggregation algorithms uses statistically linear models, including a pixel block intensity modulation algorithm [8], an emissivity modulation algorithm [9], and a high resolution urban thermal sharpener algorithm [10]. Previous work has focused on choosing other predictors such as albedo [16], percent impervious surface area [17], temperature vegetation dryness index [18], and normalized difference built-up index [19]. More recently, complex non-linear statistical algorithms with additional predictor variables have been proposed to improve performance, including least median square (LMS) regression [12], artificial neural network [20,21,22], thin plate spline interpolation [23], co-kriging method [24], wavelet transformation [25], and random forest regression [26,27]. Also, some studies have focused on spatio-temporal disaggregation that conducts data fusion between thermal imagery with low spatial and high temporal resolution and that with high spatial and low temporal resolution [28,29,30,31].

Previous work applied disaggregation methods to global applications such as evapotranspiration [20,32,33,34,35], urban heat islands [9,10,15,17,36,37], surface energy balance [38], and soil moisture estimation [39]. Practical applications of thermal data, for example, wildfires [40,41,42], volcanic activity [43,44], and land cover classification [45], however, are still limited in medium resolution sensors like Landsat-8 or advanced spaceborne thermal emission and reflection radiometer (ASTER) products. Some previous work reported that the VI-based disaggregation method of the TIR was limited in heterogeneous areas [12,22,46]; however, the limitations of previous VI-based disaggregation on fire affected areas have not been reported. The disaggregation of TIR in relatively medium-high spatial resolution, therefore, is required to expand the application of TIR products for fire incidents. Furthermore, at high temperatures up to 1500 K, temperature estimation using the TIR sensor is limited because of its radiometric limitation. In addition, surface reflectance at the VNIR bands, mostly used for the input of the disaggregation process, are distorted in wildfire situations because the presence of smoke results in low transmittance. Therefore, there is a necessity to use other indices or band imagery in the disaggregation of TIR process.

An alternative disaggregation method is required to improve previous work. There are essentially two issues with the disaggregation of TIR that require improvement. First, a new input index must be chosen such that a new disaggregation of TIR will exclude the index based on VNIR wavelengths data. After choosing a new input for the disaggregation process, there remains an issue about how to preserve the thermal variation of the disaggregated result that is similar to that of the original TIR image at coarser resolution. To deal with these issues, this study chose a SWIR band as the input of the alternative disaggregation algorithm. The SWIR band imagery has a finer resolution and can be used to estimate high temperatures in flaming combustion situations due to its radiometric feature [40]. In addition, in order to preserve spectral information in the disaggregated result, a guided image filter (GF) [47], that is, an edge-preserving smoothing operator, may be employed.

The objective of this study was to propose an efficient and improved disaggregation algorithm of TIR imagery at the scene of a wildfire using guided SWIR band imagery. The imagery was obtained from Landsat-8 operational land imager (OLI) and thermal infrared sensor (TIRS). The main parts of this paper are structured as follows. Section 2 introduces previous thermal disaggregation techniques, and the concept of a guided image filter is presented in Section 3. The proposed thermal disaggregation method is presented in Section 4, and Section 5 presents the experimental results, along with image quality assessments. Finally, a discussion of future advancements and a conclusion are provided in Section 6 and Section 7 respectively.

2. Previous Disaggregation Methods of Thermal Data

2.1. Disaggregation of Radiometric Temperature

The DisTrad method focuses on the inverse relationship between temperature and NDVI, with a 2nd order polynomial regression [6]. The disaggregation step of the algorithm is expressed by the following equations:

$${\widehat{T}}_{\mathit{LR}}=f({\mathit{NDVI}}_{\mathit{LR}})=a_0+a_1{\mathit{NDVI}}_{\mathit{LR}}+a_2{{\mathit{NDVI}}_{\mathit{LR}}}^2$$

(1)

$$\Delta {\widehat{T}}_{\mathit{LR}}=T_{\mathit{LR}}^{\mathit{obs}}-{\widehat{T}}_{\mathit{LR}}$$

(2)

$${\widehat{T}}_{\mathit{HR}}=f({\mathit{NDVI}}_{\mathit{HR}})+\Delta {\widehat{T}}_{\mathit{LR}}$$

(3)

where ${\widehat{T}}_{\mathit{LR}}$ is predicted temperature at low resolution, $T_{\mathit{LR}}^{\mathit{obs}}$ is observed temperature from satellite platform, and ${\widehat{T}}_{\mathit{HR}}$ is predicted temperature at high resolution. A linear least square regression was applied between low resolution temperature imagery ($T_{\mathit{LR}}$) and NDVI data (${\mathit{NDVI}}_{\mathit{LR}}$), and the residual term ($\Delta {\widehat{T}}_{\mathit{LR}}$) obtained by Equations (1) and (2). The low resolution regression, then, is applied to the high resolution NDVI data to estimate disaggregated temperature imagery. Finally, to reduce thermal sharpening error, the residual term from the regression at the low resolution case is applied to the estimated temperature imagery at high resolution in Equation (3).

2.2. Temperature Sharpening

The TsHARP algorithm suggested by Agam, et al. [7] has a similar disaggregation step to the DisTrad, with the difference occurring only in the sharpening basis function. In TsHARP, fractional vegetation ($f_c$) cover was selected for the input of the sharpening basis function because it is more physically correlated with temperature than NDVI. Thus, the sharpening basis function is replaced by Equation (4),

$${\widehat{T}}_{\mathit{LR}}=f({\mathit{NDVI}}_{\mathit{LR}})=a_0+a_1(1-{\left(\frac{{\mathit{NDVI}}_{\mathit{max}}-{\mathit{NDVI}}_{\mathit{LR}}}{{\mathit{NDVI}}_{\mathit{max}}-{\mathit{NDVI}}_{\mathit{min}}}\right)}^{0.625})=a_0+a_1{f}_c$$

(4)

where ${\mathit{NDVI}}_{\mathit{max}}$ and ${\mathit{NDVI}}_{\mathit{min}}$ are the maximum and minimum value of NDVI in the single scene of data. The remainder of the thermal sharpening process is equal to that of DisTrad.

2.3. Least Median Square Regression

The DisTrad and TsHARP algorithms are based on the ordinary least square (OLS) regression method. However, the OLS regression is sensitive to outliers and lacks robustness, especially in heterogeneous landscapes. Mukherjee, et al. [12] applied the regression to thermal sharpening with the NDVI in heterogeneous landscapes. Median is known as a rank statistic and is less sensitive to extreme outliers. The LMS method can be expressed as Equation (5):

$$\begin{array}{l}\mathrm{MinMedSR}\\ =\mathrm{Median}\left\{{\left(T_1-\left({a+b\times \mathit{NDVI}}_1\right)\right)}^2,{\left(T_2-\left({a+b\times \mathit{NDVI}}_2\right)\right)}^2,\dots ,{\left(T_n-\left({a+b\times \mathit{NDVI}}_n\right)\right)}^2\right\}\end{array}$$

(5)

3. Guided Image Filter

The GF is an edge-preserving smoothing operator based on a local linear ridge model using a guidance image Y as one of the inputs to filter the input image X. Originally, the GF was first suggested by He et al. [47], but there are few previous studies that applied the GF to a remote sensing, pan-sharpening process [48,49,50]. Here, we apply the GF to the disaggregation of TIR. Thus, the key features of the GF are preserving major information about the input image and reflecting variation tendencies of the guidance image [49]. The output image Z, from the GF, is assumed to be a linear transformation of the guidance image Y in a local window ${\nu }_k$, at a pixel k, as shown in Equation (6):

$${\mathit{Z}}_i{=a}_k{\mathit{Y}}_i{+b}_k{\forall }_i\in {\nu }_k$$

(6)

where $a_k$ and $b_k$ are the local linear coefficients, and ${\mathit{Y}}_i$ and ${\mathit{Z}}_i$ are the ith pixel values of the guidance image and the output. This local linear model ensures that if Y has an edge, Z also has a high spatial frequency because of the relationship that $\nabla \mathit{Z}=a\nabla \mathit{Y}$. These coefficients can be estimated by minimizing the squared difference between Z and input image X with a local ridge regularization parameter $\omega$ in Equation (7):

$$\mathrm{E}\left(a_k{,b}_k\right)={\displaystyle \sum }_{iϵ{\upsilon }_k}\left({{(a}_k{\mathit{Y}}_i+b_k-{\mathit{X}}_i)}^2+{{\omega a}_k}^2\right)$$

(7)

The regularization parameter $\omega$ determines what edges with high variance should be preserved. The local linear coefficients are given by ridge solution as shown in Equations (8) and (9):

$$a_k=\frac{\frac{1}{\left|\upsilon \right|}{{\displaystyle \sum }}_{iϵ{\upsilon }_k}{\mathit{X}}_i{\mathit{Y}}_i-{\mu }_k{\overline{\mathit{X}}}_k}{{{\sigma }_k}^2+\omega }$$

(8)

$$b_k={\overline{\mathit{X}}}_k-a_k{\mu }_k$$

(9)

where ${{\sigma }_k}^2$ and ${\mu }_k$ are the variance and mean of Y in ${\upsilon }_k$, $\left|\upsilon \right|$ is the number of pixels in ${\upsilon }_k$, and ${\overline{\mathit{X}}}_k$ is the mean of $\mathit{X}$ in ${\upsilon }_k$. The averaging strategy of overlapping windows is applied to obtain the output image Z because a pixel i is involved in all the overlapping local windows, ${\upsilon }_k$. Thus, the filtering output is written by Equation (10):

$$q_i= \frac{1}{\left|\upsilon \right|}{\displaystyle \sum }_{k|i\in {\upsilon }_k}(a_k{\mathit{Y}}_i+b_k)=\overline{a_i}{\mathit{Y}}_i+\overline{b_i}$$

(10)

where $\overline{a_i}=\frac{1}{\left|\upsilon \right|}{{\displaystyle \sum }}_{k\in {\upsilon }_i}{a}_k$ and $\overline{b_i}=\frac{1}{\left|\upsilon \right|}{{\displaystyle \sum }}_{k\in {\upsilon }_i}{b}_k$ are the average coefficients of all overlapping windows. Finally, for convenient representation of GF, this study abbreviates GF as follows in Equation (11):

$$\mathit{Z}={\mathit{GF}}_{\nu ,\omega }(\mathit{X},\mathit{Y})$$

(11)

where ${\mathit{GF}}_{\nu ,\omega }$ is a guided filter defined by filter size, $\nu$, and local ridge constant, $\omega$. The GF has useful advantages of detail enhancement. First, the GF has edge- and gradient-preserving smoothing properties. The main objective of the GF is smoothing, but the edge-preserving property prevents a loss in detail during the smoothing process, primarily near the edge. Furthermore, compared to that of a bilateral filter, also known as a smoothing operator similar to the GF, the gradient-preserving property of the GF can block the ringing or gradient-reversal artifacts near the edge. The other advantage of the GF is that the guidance image can transfer the structure itself to the input image. Therefore, the GF could be applied to SWIR imagery at a finer resolution using TIR imagery at a coarser resolution as the guidance image.

4. Proposed Thermal Disaggregation Method

There are two issues with the proposed disaggregation algorithm of the TIR data. One is the choice of a new input index, and the other is how to apply the GF in the disaggregation process of thermal products. In this study, the SWIR band was chosen as the input of the disaggregation algorithm, and the GF was adopted to estimate the geometric detail (GD) of the thermal product.

4.1. Characteristics of SWIR Wavelength

Planck’s equation can calculate the emitted electromagnetic radiance of a blackbody ($M_{\lambda ,T})$ in Equation (12):

$$M_{\lambda ,T}=\frac{{2\mathit{hc}}^2}{{\lambda }^5{(e}^{\mathit{hc}/\lambda \mathit{kT}}-1)}$$

(12)

where h is Planck’s constant, c is the speed of light, $\lambda$ is wavelength, k is Boltzmann’s constant, and T is the temperature of a blackbody in Kelvin. The wavelength of maximum emittance (${\lambda }_{\mathit{max}}$) and the total radiance for a blackbody ($M_{\mathit{tot}}$) depend on the temperature of a blackbody, and can be described by Wien’s displacement law and Stefan–Boltzmann’s law as follows in Equations (13) and (14):

$${\lambda }_{\mathit{max}}=\frac{T}{A}$$

(13)

$$M_{\mathit{tot}}={\sigma T}^4$$

(14)

where A is Wien’s constant equal to $2897.8\mathsf{\mu }\mathrm{m}\mathrm{K}$, and $\sigma$ is Stefan–Boltzmann constant equal to $5.6697\times {10}^{-8}\mathrm{W}/{\mathrm{m}}^2{\mathrm{K}}^4$. Based on both theorems, as the temperature of a blackbody increases, the wavelength of peak radiance shifts to shorter wavelengths and the total emitted radiance from the blackbody increases as shown in Figure 1, and the shorter wavelength is more sensitive to higher temperatures. In terms of emitted radiation, Figure 1 shows that the maximum radiation at a temperature of 1300 K appears at a wavelength of 2.2 μm, which means that it is highly efficient in detecting wildfires, considering land surface temperatures range from 1000 K to 1500 K under flaming combustion [40], and the SWIR-2 spectral response function of the Landsat-8 OLI includes this wavelength. However, the detected radiance by the sensor is composed of reflected solar radiation and emitted thermal radiation. According to Murphy et al. [51], a single SWIR band has been frequently saturated at a 30-m spatial resolution, and this implied that the SWIR-2 band can be affected by land cover with a high spectral reflectance in that band. Despite the solar reflective effect, the SWIR-2 band could be input into the proposed method because the longer wavelength results in a reduction of solar reflective radiance. Giglio et al. [52] reported that the longest wavelength in the ASTER SWIR band (2.4 μm) could be selected as an input considering the flaming and smoldering response in the ratio of typical fire radiance to typical land surface radiance. Therefore, the SWIR-2 band was used to obtain high spatial resolution and high frequency temperature components in flaming combustion circumstances.

4.2. Proposed Disaggregation Method

A flowchart of the proposed disaggregation method is described in Figure 2. Originally, the distributed Level 1 TIR data was already resampled to 30 m resolution by cubic convolution resampling [2]. To obtain the original resolution TIR image at 100 m resolution, an aggregation procedure was conducted on the Level 1 TIR data. However, the ratio between high- and low-resolution was rounded off to 3 because of effective image processing. The National Aeronautics and Space Administration has conducted a precise co-registration process for Landsat products [2]. Thus, the original resolution of the TIR dataset was assumed to be 90 m instead of 100 m. Also, the TIR-1 band data were selected because of the reported calibration issue. The main inputs of this method are only the TIR-1 band as brightness temperature (BT), and the SWIR-2 band as top-of-atmosphere (TOA) reflectance. After radiometric calibration, the BT imagery with 90-m resolution was upsampled at 30-m resolution, which is identical to that of SWIR imagery by cubic convolution resampling, and the upsampled BT imagery was used as guidance images for the GF. Histogram-matched SWIR imagery at 30-m resolution was selected for input images of the GF. High frequency components of BT could then be estimated as a result of GF. The difference between the TIR image and the approximation image from the GF can be considered as the GD from the TIR image. Injection gain was applied to the GD because of an energy matching issue. Finally, an injected GD of BT was merged with the original BT imagery.

4.2.1. Preprocessing

To obtain at-sensor spectral radiance, TOA reflectance, or at-sensor BT, the digital number (DN) must be converted using a radiometric calibration process [53]. First, the DN format of the image data was converted to an at-sensor spectral radiance ($L_{\lambda }$) as follows:

$$L_{\lambda }=G_{\mathit{rescale}}\times \mathit{DN}+\mathit{Offset}$$

(15)

where $G_{\mathit{rescale}}=\frac{L_{\mathit{max}}-L_{\mathit{min}}}{{\mathit{DN}}_{\mathit{max}}-{\mathit{DN}}_{\mathit{min}}}$ is a rescaling gain factor, $\mathit{Offset}=L_{\mathit{min}}-G_{\mathit{rescale}}{\mathit{DN}}_{\mathit{min}}$ is a rescaling offset factor, $L_{\mathit{max}}$ and $L_{\mathit{min}}$ are spectral at-sensor radiances that are scaled to ${\mathit{DN}}_{\mathit{max}}$ and ${\mathit{DN}}_{\mathit{min}}$, respectively, and the unit of $L_{\lambda }$ is $\mathrm{W}/({\mathrm{m}}^2\mathrm{sr}\mathsf{\mu }\mathrm{m})$. Based on the calculated at-sensor spectral radiance, TOA reflectance (${\rho }_{\lambda }$) and at-sensor BT ($\mathit{BT}$) can be obtained by employing Equations (16) and (17):

$${\rho }_{\lambda }=\frac{\pi ·L_{\lambda }·d^2}{{\mathrm{ESUN}}_{\lambda }·{\cos \theta }_s}$$

(16)

$$BT=\frac{K_2}{\mathit{ln}\left(\frac{K_1}{L_{\lambda }}+1\right)}$$

(17)

where $d$ is earth-sun distance in astronomical units, ${\mathrm{ESUN}}_{\lambda }$ is solar irradiance, ${\theta }_s$ is solar zenith angle, and $K_1$, $K_2$ are calibration constants. In this study, the DN of VNIR and SWIR imagery were converted to TOA reflectance by using the fast line-of-sight atmospheric analysis of spectral hypercubes (FLAASH) model [54] which has been incorporated into the environment for visualizing images (ENVI) software. This allows for the DN of VNIR and SWIR images to be converted to TOA reflectance and that of TIR to be converted to BT. Converted TOA reflectance and BT were then used as an input of both the proposed and comparative thermal sharpening algorithms. The use of TOA reflectance reduced the effect of reflective solar radiance following the Earth-sun-satellite geometry [51].

4.2.2. Disaggregated Geometric Detail Using GF

After the radiometric calibration, the SWIR reflectance image (${\rho }_{\mathit{SWIR}}$) at 30-m resolution was histogram aligned with BT data (BT) from the TIR sensor at 90-m resolution to match the mean and standard deviation of the BT. Because of the inherent features of SWIR wavelengths, especially the SWIR-2 band in the Landsat-8 OLI sensor, the matched SWIR reflectance (${\rho }_{\mathit{SWIR}}^{\mathit{BT}}$) was used to obtain high spatial resolution and a high frequency temperature component in flaming combustion circumstances instead of NDVI. In order to estimate the disaggregated GD (GD) of the BT image, the matched SWIR reflectance was selected as an input image of the GF, and the upsampled BT imagery was used as a guidance image of the GF. Therefore, the disaggregated GD of the BT was estimated by the GF because the guided SWIR (${\mathit{GF}}_{\nu ,\omega }{(\rho }_{\mathit{SWIR}}^{\mathit{BT}},\mathit{BT})$) can preserve high frequency temperature components, as in flaming combustion conditions, with the variance tendency of the BT. The disaggregated GD of the BT was obtained by the difference of the histogram-matched SWIR and the guided SWIR in Equation (18):

$$\mathit{GD}={\rho }_{\mathit{SWIR}}^{\mathit{BT}}-{\mathit{GF}}_{\nu ,\omega }{(\rho }_{\mathit{SWIR}}^{\mathit{BT}},\widetilde{\mathit{BT}})$$

(18)

where ${\rho }_{\mathit{SWIR}}^{\mathit{BT}}$ is histogram-matched SWIR reflectance imagery by the BT data, and $\widetilde{\mathit{BT}}$ is the upsampled BT. Thus, high frequency temperature components from SWIR images and structures of TIR images were preserved in the guided image, and the guided image could be used to obtain the disaggregated GD. Although the method of the disaggregation of thermal products influences the result of that method, the excessive or deficient window size also results in redundant information and produces artifacts in the resultant image [49]. In this study, the regularization parameter $\omega$ was set to 1, and performance with different window size ${\upsilon }_k$ was evaluated for the proposed method with quantitative image quality indices.

4.2.3. Injection Gain-Based on Data Statistics

The appropriate injection gain is an important factor for addition of GD information to the TIR image. The extraction of high frequency components and their effective injection are, generally, recognized as a single algorithm in an image fusion field [55]. If the injected spatial detail is not enough or is excessive, the spatial quality of the disaggregation image is not sufficient, or redundant information can lead to spectral distortion. In other words, in this study, the injection gain was considered as a process to reduce solar-reflected radiance in the SWIR-2 because the GD was obtained from the input of that band. Thus, this study strived to have the injection gain properly reflect the GD. Although the GD had a temperature component because the input of the GD was histogram-matched with the BT imagery, the distribution of the GD is highly skewed relative to the original BT imagery due to a high sensitivity of SWIR at the temperature of combustion. In order to achieve the proper energy matching between the GD and the BT imagery, injection gain was utilized in this study. The proposed injection gain is based on the ratio of the multiplication of the range and the skewness between the disaggregated GD and the original BT data as follows:

$$\mathsf{\Omega }=\frac{{\mathit{range}}_{\mathit{BT}}·{\mathit{skewness}}_{\mathit{BT}}}{{\mathit{range}}_{\mathit{GD}}·{\mathit{skewness}}_{\mathit{GD}}}$$

(19)

where $\mathsf{\Omega }$ is the injection gain proposed in this study, ${\mathit{range}}_{\mathit{BT}}$ and ${\mathit{range}}_{\mathit{GD}}$ are the range of the BT and the disaggregated GD, and ${\mathit{skewness}}_{\mathit{BT}}$ and ${\mathit{skewness}}_{\mathit{GD}}$ are the skewness of the BT and the disaggregated GD. Finally, considering the proposed injection gain, the disaggregation of the BT imagery can be described by Equation (20):

$$\widehat{BT}=\widetilde{BT}+\mathsf{\Omega }\odot \mathit{GD}$$

(20)

where $\widehat{BT}$ is the disaggregated BT imagery at a 30 m resolution, and $\odot$ denotes a pixelwise multiplication operator.

4.3. Image Quality Assessment

To evaluate the quantitative image quality, five widely used quality indices were selected in this study. First, the root mean square error (RMSE) and the mean absolute error (MAE) are used to fundamentally evaluate the result of the image fusion process. The RMSE and MAE represent the difference and the bias between the fused image F and the reference image R. Both can be expressed by the following Equations (21) and (22):

$$\mathrm{RMSE}=\sqrt{\frac{1}{M·N}{\displaystyle \sum }_{m=1}^M{\displaystyle \sum }_{n=1}^N{(\mathit{F}\left(m,n\right)-\mathit{R}\left(m,n\right))}^2}$$

(21)

$$\mathrm{MAE}=\frac{1}{M·N}{\displaystyle \sum }_{m=1}^M{\displaystyle \sum }_{n=1}^N(\mathit{F}\left(m,n\right)-\mathit{R}\left(m,n\right))$$

(22)

where F(m,n) and R(m,n) represent the pixel values of the fused image, F, and the reference image R, at the mth column and nth row, and M × N is the size of the reference image. Lower values for RMSE and MAE indicate a lower error for the fused image. Second, the correlation coefficient (CC) is also used to assess the fused images. The CC of a fused image and its reference image reflects the similarity of spectral information. The formulation is Equation (23) as follows:

$$\mathrm{CC}=\frac{{{\displaystyle \sum }}_{m=1}^M{{\displaystyle \sum }}_{n=1}^N\left(\mathit{F}\left(m,n\right){-\mu }_F\left)\right(\mathit{R}\left(m,n\right) - {\mu }_R\right)}{\sqrt{{{\displaystyle \sum }}_{m=1}^M{{\displaystyle \sum }}_{n=1}^N{\left(\mathit{F}\left(m,n\right) - {\mu }_F\right)}^2{{\displaystyle \sum }}_{m=1}^M{{\displaystyle \sum }}_{n=1}^N{\left(\mathit{R}\left(m,n\right) - {\mu }_R\right)}^2}}$$

(23)

where ${\mu }_F$ and ${\mu }_R$ are the mean of the fused image and the reference image, respectively. If the two images are correlated, the CC is close to 1, which implies that the spectral features of the original image were preserved well during the disaggregation process. Also, the image quality Q index, also known as the universal image quality index (UIQI), reflects the structural distortion degree [56]. UIQI is defined in Equation (24) as follows:

$$\mathrm{UIQI}=\frac{{\sigma }_{\mathit{FR}}}{{\sigma }_F{\sigma }_R}·\frac{2\left|{\mu }_F\right|\left|{\mu }_R\right|}{{\left|{\mu }_F\right|}^2+{\left|{\mu }_R\right|}^2}·\frac{{2\sigma }_F{\sigma }_R}{{{\sigma }_F}^2{{+\sigma }_R}^2}$$

(24)

where ${\sigma }_{\mathit{FR}}$ is the covariance matrix between the fused image and the reference image, and ${\sigma }_F$ and ${\sigma }_R$ are the standard deviations of the fused image and the reference image. The UIQI was applied to the disaggregated results with an 8 × 8 distinct window. If its value is close to 1, the quality of the fused image is adequate. Finally, erreur relative globale adimensionelle de synthèse (ERGAS) is useful to evaluate the overall quality of the disaggregated image [3]. It represents the degree to which the disaggregated image and the reference image are different. The ERGAS is defined in Equation (25) as follows:

$$\mathrm{ERGAS}=\frac{100}{\beta }\sqrt{\frac{1}{N}{\displaystyle \sum }_{i=1}^N{\left(\frac{{\mathit{RMSE}}_i}{{\mu }_i}\right)}^2}$$

(25)

where $\beta$ is the scale ratio between the pixel sizes of the SWIR image and the BT image, ${\mathit{RMSE}}_i$ is the RMSE between the ith fused band and the reference band, and ${\mu }_i$ is the mean of the ith reference band. Thus, a small value of ERGAS means a small spectral distortion is present in the disaggregated image.

There are two validation methods that are commonly used to assess the quality of the fused image. The first is the synthesis property, and the other is the consistency property as shown in Figure 3. First, the synthesis property considers that the disaggregated image should be as close to identical to the acquired input image as possible. However, at the full-resolution scale, assessment of the synthesis property is not possible, thus the fusion images, SWIR and TIR, were degraded by the scale ratio between the pixel sizes of the SWIR image and the BT image. The degraded images were used as an input of the disaggregation method, and the degraded fused image was compared to the original observed BT image as a reference. Second, the consistency property reflects the degree to which the spatially degraded fused image is identical to the original observed BT image. Cubic convolution interpolation was applied for image degradation of the synthesis and consistency properties. To validate the evaluation process, both synthesis and the consistency properties were applied to our experimental assessments.

5. Experimental Results

5.1. Study Area

On 6 May 2017, a series of wildfires that destroyed more than thirty homes and caused the evacuation of 2500 residents occurred simultaneously in Samcheok, Gangneung and Sangju, cities in South Korea. Together, these fires burned nearly 150 ha. Especially dry climate conditions in spring promoted the rapid growth of the wildfires. The study area is located in Samcheok city, Gangwon province, one of the areas affected by the spring 2017 fires (Figure 4). For the proposed thermal disaggregation and comparative processes, freely accessible Landsat-8 Level-1 OLI and TIRS datasets, with a spatial resolution of 30 m in the VNIR and SWIR bands and 30 m in the TIR bands (originally, 100 m), were acquired from the USGS Earth Explorer website [57]. The datasets were acquired at 01:58:02 GMT (10:58:02 local time), May 8th, 2017 from the 114 path and 34 row in the world reference system-2 (WRS-2). The acquisition datasets were clipped including the wildfires and marginal cloud-free areas, and the size of the TIR images was 267 × 543 pixels, and that of the MS images was 801 × 1629 pixels. Considering the two-day gap between the initiated wildfire and the image acquisition date, the datasets are expected to include a variety of wildfire components, for example, flaming combustion, smoldering combustion, and burned areas. In particular, flaming combustion, a core of wildfires, is remarkably shown in the acquired imagery (refer to Figure 4).

5.2. Analysis of the Disaggregation Methods

Experiments were conducted to evaluate the performance of the proposed disaggregation of TIR imagery using Landsat-8 datasets. The proposed algorithm was compared with DisTrad [6], TsHARP [7], and the LMS method [12]. Evaluation of the disaggregation methods focused on the synthesis and consistency properties [55]. The reflectance of the SWIR-2 band from the Landsat-8 OLI sensor was selected as the input of the proposed algorithm, and the GF filter constrained the input imagery by using the guidance image, and original TIR data at a coarser resolution. The results of the proposed disaggregation algorithm and the comparison disaggregation methods were considered, and the disaggregated image quality assessment was then conducted with parameter configuration (window size) of the GF.

5.2.1. Visual Analysis

The proposed disaggregation method was compared with the DisTrad, TsHARP, and the LMS methods. Figure 5 shows the result of the disaggregated BT image obtained using the proposed SWIR and the GF-based algorithm with the color composite image, the original BT image at a resolution of 90 m, and the distributed BT image from Earth Explorer at a 30 m resolution. To qualitatively evaluate the results, three subsets of the disaggregated images, A, B, and C in Figure 5 were obtained, and comparison with VI-based methods including DisTrad, TsHARP, and LMS was conducted directly in the selected subset areas as shown in Figure 6, Figure 7 and Figure 8. There were two notable features of the results from the proposed method in terms of qualitative analysis. First, the proposed algorithm performed well on wildfire areas, spatially preserving high frequency temperature components (Figure 6 and Figure 7). The comparison algorithms had lower BT values in the disaggregated image on burning areas, but the proposed algorithm produced higher BT values on the same burning areas. Also, removal of ringing artifacts was demonstrated in the disaggregated image of our proposed algorithm. Figure 8 shows the subset of the disaggregated images near the lake. The ringing artefacts near the lake edge were eliminated by the proposed algorithm but appear in the results of other algorithms (Figure 8). As shown in the visual evaluation, this experiment confirmed that the proposed algorithm can provide visually satisfactory results with regard to the original BT image.

5.2.2. Quantitative Analysis

Table 1. Quantitative evaluation of the proposed method depending on different window size. The optimal values in terms of quality index are highlighted in bold.

	Window Size
	3 × 3	5 × 5	7 × 7	9 × 9	11 × 11	13 × 13	15 × 15	17 × 17	19 × 19
RMSE	0.073	0.072	0.073	0.074	0.074	0.075	0.076	0.076	0.077
CC	0.9996	0.9996	0.9996	0.9996	0.9996	0.9996	0.9995	0.9995	0.9995
MAE	0.051	0.049	0.049	0.049	0.049	0.050	0.050	0.051	0.051
UIQI	0.9984	0.9985	0.9985	0.9985	0.9984	0.9984	0.9984	0.9983	0.9983
ERGAS	0.025	0.025	0.025	0.025	0.025	0.026	0.026	0.026	0.026

presents the results of quantitative analyses depending on different window size. As shown in Table 1, the best performance of the proposed method was obtained from window size 5 × 5 of the GF. The units of RMSE and MAE are Kelvin in this study. The results of the disaggregated image quality assessment based on the synthesis property are presented in

Table 2. Results of quantitative image quality assessment corresponding to each disaggregation algorithm. The optimal values in terms of quality index are highlighted in bold.

	Synthesis Property					Consistency Property
	RMSE	MAE	CC	UIQI	ERGAS	RMSE	MAE	CC	UIQI	ERGAS
DisTrad	0.796	0.513	0.9515	0.8414	0.271	0.175	0.100	0.9979	0.9940	0.060
TsHARP	0.609	0.463	0.9700	0.8261	0.208	0.100	0.077	0.9992	0.9953	0.034
LMS	0.607	0.460	0.9703	0.8282	0.208	0.099	0.075	0.9992	0.9955	0.034
GF	0.472	0.331	0.9822	0.9294	0.161	0.072	0.049	0.9996	0.9985	0.025

.

The disaggregated image obtained by the DisTrad algorithm indicates that all quality indices had poor results, but the TsHARP algorithm, a refined version of the DisTrad, shows better results in comparison with the DisTrad. The disaggregation result of the LMS algorithm was similar to that of the TsHARP algorithm. Finally, the results of the proposed algorithm, indicated by GF, have the best values of RMSE, MAE, CC, UIQI and ERGAS compared to other algorithms. The thermal variance of the original BT data was therefore preserved in the disaggregated result for the proposed GF algorithm. For the consistency property, the results of quantitative analysis also show similar trends compared to the synthesis property, and here the proposed GF algorithm provides the best quality indices results for all consistency properties (Table 2).

6. Discussion

In this study, disaggregation of BT imagery of Landsat-8 from 90-m resolution to 30-m resolution on wildfire areas was conducted by using the guided SWIR-2 band of Landsat-8 OLI. The experimental results showed that the comparison methods were limited on wildfire areas, but the proposed algorithm yielded better disaggregated BT images in terms of qualitative and quantitative analysis. In particular, the qualitative evaluation of the experimental results demonstrated that the high temperature component from SWIR was preserved in the guided filtering process on the burning areas, and the SWIR band as the input of disaggregation process supplemented the limitation of TIR sensor of restricted detection of high temperature in active fires. Also, most previous disaggregation studies improved spatial detail from the MS band products, but did not consider the thermal variation of the TIR data. This occurrence of spectral distortion in the disaggregation process results because of an injection effect of redundant detail. In contrast, the thermal variation of the original BT image was preserved simultaneously with improved spatial resolution of the BT by the GF, using the original BT data as a guidance image as described in the quantitative analysis discussion (Section 5.2.2).

In the study area, there are two main wildfire components, flaming and smoldering. Despite the mixed-pixel problem within the wildfire component and unburned fuel component, Yu et al. [58] reported that the critical fraction of recognition was 1/500 of the 30-m resolution cell on an experimental hotspot target (955 K) in the Landsat SWIR band, and Giglio et al. [52] suggested that the TIR band (8.3 μm) of the ASTER data would be the best choice for detecting smoldering considering the solar reflective effect. As a result, for a flaming combustion region, the disaggregated BT was relatively higher in the margin area in the proposed disaggregation method due to the input of the SWIR-2 band, and the smoldering effect was also reflected by the tendency of the TIR band. However, temperature of the combustion region is not equal to the disaggregated BT due to a land surface temperature range from 1000 K to 1500 K in flaming combustion conditions [40]. Also, conversion the BT into land surface temperature was not conducted in this work. Although these limitations existed in this study, the application of the proposed disaggregation algorithm was proved valuable for middle-high spatial resolution imagery using Landsat-8 products, considering the limited Landsat-8 data availability on the active fire areas.

Despite applications of the TOA reflectance as an input and injection gain, the disaggregation results could hardly be free from reflective solar radiation for daytime analysis. Manmade land cover (e.g., concrete, rooftops, and highly reflective paint) and arid soil area especially have high reflectance in the SWIR-2 band. The proposed algorithm additionally used the TIR imagery, which is insensitive to reflected solar radiation [59], during the process of GF filtering. In this case, the TIR imagery was used as a guidance image and compensated for the high frequency components from the high solar reflectance. Also, the effect of reflective solar radiation depends on the characteristics of the experimental area. The study area is a typically rural and mountainous area with a low population density. In terms of disaggregation, manmade land cover could be a high frequency temperature component during the daytime because the temperature of manmade structures is higher than that of vegetation, and the high solar reflectance reflected the tendency of the BT in the experimental areas. Only for a high solar reflectance near the flaming case could the high temperature term be extended near the flaming area. However, the thermally emitted radiation must be distinct from the solar-reflected radiation in the case of an arid or urban area. To avoid a saturation problem, an alternative daytime hot spot detection algorithm was suggested using separate wavebands including NIR, SWIR-1, and SWIR-2 of the Landsat datasets [51]. Therefore, additive band datasets, the NIR and the shorter SWIR band, will be considered in a future work to reduce the effect of solar reflectance and potential saturation in the SWIR-2 band.

This study expected to apply the proposed algorithm to improve the spatial resolution without spectral distortion in the production of Landsat thermal images Application of the proposed method can be extended to geostationary products that have low spatial and high temporal resolutions, for example, moderate resolution imaging spectroradiometer (MODIS) and visual infrared imaging radiometer suite (VIIRS) which cover the SWIR and TIR bands with diurnal temporal resolution. Consequently, data fusion with the advantages of high temporal geostationary products and high spatial resolution, termed spatial-temporal disaggregation, can be conducted using the proposed method based on the GF. Future work will be focused on the application of the proposed disaggregation algorithm to improve the spatial-temporal resolution of the TIR data.

7. Conclusions

In this study, the efficient and improved disaggregation algorithm of TIR imagery using guided SWIR imagery was proposed for wildfire conditions. The SWIR imagery radiometrically satisfied the preservation of a high temperature component in flaming combustion, and the GF prevented the spectral distortion of the original TIR data and simultaneously preserved the high frequency component of SWIR image in the disaggregated result. Landsat-8 OLI and TIRS images, captured at the wildfire site, were used to evaluate the proposed algorithm with the existing disaggregation algorithms based on VI originated from VNIR. The experimental results demonstrated that the proposed algorithm disaggregated the BT data from 90-m resolution to 30-m resolution and satisfied image quality assessments in terms of visual and quantitative evaluation. In the quantitative evaluation, RMSE, MAE, CC, ERGAS, and UIQI had better results (Table 2) in the proposed algorithm than those of the other comparison methods. Also, visual analysis (Figure 6 and Figure 7) demonstrated that the proposed algorithm conducted the disaggregation of the TIR imagery reflecting the temperature of wildfire areas. The experimental results show that the proposed algorithm was successfully applied to Landsat-8 OLI and TIRS datasets from wildfire areas. In future work, the proposed method will be developed considering the effect of SWIR saturation, and also applied to spatio-temporal disaggregation with a geostationary image and Landsat datasets.