Improving remote sensing of extreme events with machine learning: land surface temperature retrievals from IASI observations

The IASI instrument is a Michelson interferometer measuring infrared radiation placed onboard the European meteorological satellites MetOp. The MetOp series is operated by the European Organisation for the Exploitation of Meteorological Satellites (EUMETSTAT). IASI offers a high resolution infrared spectrum between 675 and 2760 cm⁻¹ over 8461 channels. This spectrum includes several window channels for the retrieval of surface properties, from which we regularly extract channels between 770–980 cm⁻¹, 1080–1150 cm⁻¹ and 1210–1650 cm⁻¹ (with a total of 563 channels). This extraction is a good compromise that includes most of the information contained in these window channels while keeping some redundancy to help reduce the instrument noise. We then apply a principal component analysis (PCA), an effective method to reduce dimensionality [17, 27–29] prior to using the ML/DL algorithms. The information is compressed using a PCA into three principal components (PCs) representing 99.30%, 0.54% and 0.08% of the total variance (for a total of 99.92%). These PCs contain most of the information about surface properties from the IASI spectrum. Such a high percentage for the primary PC can be explained by the high spectral correlation between the extracted channels. We extract an image each time a MetOp A polar orbit partially covers the chosen domain (approximately twice a day), which results in 15 711 'images' over the period. Only 8242 images of the domain over the ten-year period are kept after screening out the images only covering a small part of the domain or those with a large cloud fraction. The acceptance criteria is that there should be at least 30 clear pixels in the image, and a Euclidean distance of at least 20 pixels between a pair of clear pixels. Naturally, the extracted images only contain clear-sky pixels, since no information on the surface can be obtained for cloudy pixels. An analysis of the extracted database was performed to ensure that both summer and winter months were fairly represented.

The ECMWF ERA5 reanalysis [25] is a re-processing of an NWP model constrained with most of the available observations. It is used to estimate the state of the surface and atmosphere with an hourly time-step. The LST product over the domain is extracted from this reanalysis, and the observed clear-sky IASI TBs are collocated with the closest hour of the LST product.

A climatology of surface emissivities derived from TELSEI [18, a Tool to Estimate Land Surface Emissivity in the Infrared] is used. This tool gives monthly global estimates of infrared emissivity with IASI spectral resolution. To stay coherent with the extracted IASI TB channels, three channels are kept for the emissivity database. As in [15], we choose here the 850, 900 and 1100 cm⁻¹ wave-numbers to be representative of the overall spectrum.

MLPs are a simple and very common type of NNs, in which neurons of one layer have direct connections to neurons in the next layer. If enough samples are available, and the architecture contains enough neurons, any continuous relationship can be learnt by an MLP model [35–37]. An MLP tries to approximate a function $f(x) = y$ as $\widehat{f}(x) = \widehat{y}$, where x are the inputs to the MLP and $\widehat{y}$ the MLPs approximation of the targets y. The quality of $\widehat{y}$ strongly depends on the information available in the inputs x. In the remote sensing community, the inputs to MLP models are generally from single pixels only. This is both because there is a common consensus that there is little or no information that can help the retrieval among the neighboring pixels, and also because the nature of the data (i.e. shifting orbits) does not facilitate the use of image-processing approaches.

We nonetheless propose to test CNNs, one of the most popular image-processing models, to exploit the information in neighboring pixels and spatial patterns to improve the quality of retrievals. CNNs are deep NN models where each convolutional layer is a set of filters passed over the height and width of the input volume (typically, lat × lon × PCs). After training of the NN is complete, filters activate when they detect a certain spatial pattern. These approaches require to use images in input to the model. For our case study, it is therefore necessary to transform the original IASI data (i.e. orbits) into images. Here, this is done by focusing on a specific domain and capturing an image every time an orbit passes over the chosen area. Other solutions for IASI data are considered in [16]. As previously mentioned, only clear-sky pixels can be used to retrieve LST. However, CNNs require complete images, both for the training and operational prediction phases. This is done by inter- and extrapolating the images. The remaining missing pixels are considered larger holes. These holes are filled using an algorithm relying on the empirical orthogonal function (EOF) representation of images [27, 28, 38]. In this approach, an image X can be rewritten as $X = \lambda_1 \cdot F_1 + \lambda_2 \cdot F_2 + \ldots + \lambda_m \cdot F_m$. After calculating the EOFs ($F_1,{\ldots}F_m$) using the covariance matrix of our database (computed using only clear pixels) the missing data are filled with a first guess: in this case we use $\lambda_1 \cdot F_1$ with λ₁ that minimizes the root mean square error over the available pixels. To finish, an optimization procedure is used to project the now complete image onto the EOF basis to obtain $\lambda_1 \ldots \lambda_m$ and replacing the missing data with $\lambda_1 \cdot F_1 + \lambda_2 \cdot F_2 + \ldots + \lambda_m \cdot F_m$, until convergence of the λ coefficients. This protocol is described more in depth in [15].

We introduce the dampening effect, defined by the range of values in the retrieved LST dataset being smaller than the range of values in the original (ERA5) LST dataset. Similarly, we introduce the inflating effect as the opposite; the range of values in the retrieved LST dataset being larger than the range of values in the original dataset.

Mathematically, let us define the set of samples of extreme low and high values of LST at the pixel level (pixel i) from the 10th lower ($Q_{0.1}^\mathrm{ERA5}(i)$) and higher ($Q_{0.9}^\mathrm{ERA5}(i)$) percentiles:

$$\begin{equation} S_\mathrm{L}(i)= \{y^t\ |\ y^t_\mathrm{ERA5}(i) \leqslant Q_{0.1}^\mathrm{ERA5}(i) \} \textrm{ and } \end{equation} \tag{ 1 }$$

$$\begin{equation} S_\mathrm{H}(i) = \{y^t\ | \ y^t_\mathrm{ERA5}(i) \geqslant Q_{0.9}^\mathrm{ERA5}(i) \}, \end{equation} \tag{ 2 }$$

where t spans the time dimension. This definition is valid both for the original (ERA5) and retrieved dataset. We then define an estimate of the range over a pixel i, for a particular dataset, as:

$$\begin{equation} R_i =\left\{\overline{S_\mathrm{H}(i)} - \overline{S_\mathrm{L}(i)}\right\}. \end{equation} \tag{ 3 }$$

The range R_i can thus be understood as the average LST value of samples in $S_\mathrm{H}(i)$ minus the average LST value of samples in $S_\mathrm{L}(i)$. Dampening is the process by which the retrieved range, computed from the $\widehat{y}$ values, is smaller than the real range, computed from the $y_\mathrm{ERA5}$ values:

$$\begin{equation} R_{i,\mathrm{RET}} - R_{i,\mathrm{ERA5}} < 0. \end{equation} \tag{ 4 }$$

Because the retrieval errors are larger for these extreme values, it is often argued that the problem is due to a lack of extreme conditions in the training dataset. However, enriching the training dataset of the NN with additional extreme conditions does not solve the problem because sampling is not the real issue.

Figure 2(a) aims to illustrate what happens. We consider a NN with only one input (TB PC1) and one output (LST). Let us consider a cloud of points (in red) representing the actual true input-to-output relationship (i.e. the TB-to-LST relationship). As previously described, an NN is an algorithm that tries to approximate a complex function $f(\hbox{TB}) = \hbox{LST}$. The NN can be seen as a injective function with only one output for each input. This function is shown in blue in the figure. Since the relationship is monovariate, the NN approximates the cloud of points for every little horizontal bin along the x-axis with an average value of the target LST (the black squares on the figure). It can be seen that the NN approximation (in green) runs through these points. The extreme values (both smallest and largest) in the shaded areas will thus not be present in the retrieval range. Adding samples to represent these extreme values will not reduce the vertical dispersion. In fact, no matter the statistical method used (linear regression, random forests, NN, etc), this vertical dispersion cannot be reduced. The only possible solution would be to introduce more PCs into the model with the expectation that the additional information will help explain and predict the existing vertical dispersion in all horizontal bins.

Zoom In Zoom Out Reset image size

When working on a spatial domain like it is the case in this study, a form of data pooling is commonly done. Data pooling refers to combining multiple pixels of the domain to train a model, over which the TB-to-LST relationships differ. This data pooling creates an overall dampening on the global range. However, at a local scale, the LST ranges can in fact be wider than what they should be. This happens because the obtained NN performs a compromise among all the pixels. The expectation when using a global NN is that the retrieval performs well on each location (i.e. over the Alps, over coastal areas $\ldots$). This is an ambitious task, all the more so when working at the global scale. This would be possible if we had all the necessary information to correctly represent the radiative transfer (RT) (forward or inverse) all around the world, whatever the environment type or local conditions, but it is difficult for the Earth surfaces where the RT is still not entirely satisfactory [39]. Some surface parameters describing the soil properties, or the state of vegetation are not available at the global scale so only a simplified relationship is attainable. A global NN would in fact find the best compromise, using a simplified relationship between the satellite observations and the variable(s) to retrieve. Due to this, we may obtain an inflating effect:

$$\begin{equation} R_{i,\mathrm{RET}} - R_{i,\mathrm{ERA}} > 0. \end{equation} \tag{ 5 }$$

This inflation is contradictory with the nature of regression models. This is why there are fewer pixels with inflating than dampening. Figure 2(b) shows the dampening and inflating effects by plotting the probability density function (PDF) of the difference in ranges ($\Delta_R = R_\mathrm{RET}-R_\mathrm{ERA5}$) for all pixels. There is a net dampening when $\Delta_R \lt 0$, and a net inflating when $\Delta_R \gt 0$.

Both the dampening and inflating effects have negative consequences on the retrieval (i.e. they contribute to the error of the model), especially on extreme events. Whether one wishes to focus on minimizing one or the other depends on the particular application. For instance one may choose to rather focus on reducing the dampening error to not miss an extreme event or the inflating effects to avoid false alarms. The goal in this study is to improve the estimation of extreme cases by limiting both the dampening and inflating effects as much as possible.

Localization, defined as a way to help the NN adjust its behavior to local conditions [15], is a powerful tool to improve retrievals [16]. Localization of NNs is a new idea in the remote sensing community. Indeed, the majority of NN models used nowadays in remote sensing are global and standard CNN models use weight-sharing. However, it was shown in [15, 16] that using localization techniques improves the NN results by adjusting its behavior to specific local conditions. This is very important for operational retrievals from satellite observations. Localization limits the need for data pooling and will therefore limit the inflating effect. In addition, localization improves the retrieval across all horizontal bins, and should therefore limit the dampening effect too. Boucher et al [15, 16] also discuss the advantages and disadvantages of localization.

Several localization approaches exist, both for MLP models and image-processing techniques like CNNs. A possibility is to add more information into the NN inputs. These auxiliary variables help the NN to understand better the physical conditions on each pixel, by better constraining the RT on this location. It is therefore able to better adapt its behavior to varying environments. Another form of auxiliary variables was introduced in [40], in which visible bands and thermal images are combined, to propagate high-resolution information for the retrieval of SST. In this study, we solely use surface emissivities as a localization variable [18, see section 2.1.3]. This approach is possible for both MLP and CNN models.

Alternatively, it is possible to use independent MLP models for every pixel of the domain. This way, the models do not have to make any compromise since they focus only on the input-to-output relationship over one particular pixel. By localizing using pixel-based NNs, less information is required to fully describe the RT on each location: some auxiliary variables might be omitted, in particular if they are constant on each pixel.

Another localization consists in a modified version of CNNs: the Localized-CNN, in which the convolution layer is thus replaced with a locally connected layer, in which specialized spatial filters are derived for every pixel of the domain instead of using the same filters everywhere (associated to weight sharing). Chen et al [41] introduce this layer (that we will call a local convolution) as an intermediate solution between a fully connected layer and a classic convolutional layer. Using such a layer for remote sensing applications was first introduced in [15, 16] and was proven to be a successful form of localization to improve retrieval statistics. In the following, several retrieval methods with different localization strategies will be compared:

• A Unique MLP model with TB PCs and emissivities in the inputs. This model is composed of one hidden layer with 100 neurons followed by a linear activation function. The Levenberg–Marquardt optimization algorithm is used for training.
• Independent MLP models for each pixel of the domain. A separate model is trained over each pixel of the domain, composed of one hidden layer with five neurons followed by a linear activation function. The inputs to this model are the TB PCs. The Adam optimization algorithm [42] is used for the training.
• A Localized-CNN model with TB PCs and emissivities. The architecture of this model is composed of a single localized-convolutional layer with a single $5\times5$ kernel. Experiments demonstrate that one convolutional layer with one kernel is enough to capture the spatial dependencies present in the full image. This is due to the highly localized property of the layer. This results in a unique feature map. The Adam optimization algorithm [42] is used for the training.

Training hyper-parameters for each model are shown in table 1. Please note that different optimization algorithms are used because in the case of the unique MLP, the problem is more complex: a lot of different local conditions are pooled into one big database and optimizing the MLP model with the Adam algorithm did not lead to satisfactory results. The Levenberg–Marquardt algorithm is more efficient in this case. However, this algorithm was not available to train CNNs in Python, therefore we opted for the well-known Adam algorithm in this case. Furthermore, all models are trained with separate training and validation sets and early stopping was used to ensure that no overfitting occurs. The training time of the Unique MLP model is 5 min. Each independent MLP model takes 1–2 s to train, depending on the pixel location. Lastly, the localized-CNN model takes 8 min to train.