Improving the Thermal Infrared Monitoring of Volcanoes: A Deep Learning Approach for Intermittent Image Series

Nine active volcanoes were targeted in this project: Mount Erebus (Antarctica), Erta Ale (Ethiopia), Mount Etna (Italy), Kīlauea (Hawaii), Masaya (Nicaragua), Nyamulagira (Democratic Republic of Congo), Nyiragongo (Democratic Republic of Congo), Pacaya (Guatemala), and the distant Pu’u ’Ō’ō cone of Kīlauea (Hawai’i, US). These volcanoes were chosen for their active nature and/or because they have displayed lava lakes (a superficial and persistent thermal phenomenon) for at least some percentage of the data record. Among these, the Hawaiian volcanoes and Mount Etna are among the most studied and monitored in the world; Erebus, Masaya, and Pacaya have some degree of ground monitoring; Nyamulagira and Nyiragongo have decommissioned ground monitoring efforts due to civil unrest in the region; and, Erta Ale does not have regular ground monitoring [2].

Among thermal monitoring options, nighttime imagery is preferred in the field of volcanology because it is generally less cloudy and decreases the effects of topography and solar heating, “allowing better delineation of low-temperature geothermal anomalies” [43, p. 18]. In general, many platforms do not make nighttime thermal imagery easily accessible (e.g., Landsat has special request targets [44]); however, the $>20$-year NASA Terra platform which hosts the MODIS and ASTER instruments are an exception to this. MODIS is very popular because it provides a very high temporal sampling rate (four images per day [45]) which can aid in the rapid detection of changes [29]. ASTER, on the other hand, provides a 10-fold increase in spatial resolution which allows the less frequent detection of subtler changes months in advance [46, 47, 5]. Nighttime thermal data from ASTER (L1T) is available for bulk download in the form of 5 bands of radiance values (8-12$\mu$m) from the online DAAC2Disk Utility, courtesy of the NASA Land Processes Distributed Active Archive Center (LP DAAC), USGS Earth Resources Observation and Science (EROS) Center, Sioux Falls, South Dakota, https://lpdaac.usgs.gov/tools/daac2diskscripts/. The level-2 Surface Kinetic Temperature product (AST_08) [48] which is more thoroughly corrected for the atmosphere and with radiance converted to temperature is available for order from the NASA LP DAAC via NASA’s Earthdata Search (https://search.earthdata.nasa.gov/search?q=ast_08).

First, all L1T data for each volcano between the dates December 18, 1999, and May 1, 2020, were downloaded. Initial inspection revealed a poor correlation between viable scenes and cloud cover metadata, so all scenes were manually inspected for viability in the data set. Upon manual inspection, each scene was denoted viable, uncertain, or nonviable (Figure 1). Due to the large number of nonviable scenes, uncertain scenes were ultimately included in our data set; these were scenes with partial cloud cover or scenes with large extents of cloud cover, through which, thermal features were still observable. The final set of viable and uncertain scenes were then ordered as level-2 AST_08 data for the remainder of the project. The distribution of time between these scenes is presented in Figure 2.

The AST_08 data contained missing values representing areas outside the sensor’s field of view as well as pixels amid very volcanically active areas (see Figure 3, top row) – with sources referencing the latter as recovery pixels [49]. In the level-2 data, the latter appears as counter-intuitive dark spots within a bright hotspot, while the level-1 data does not contain these missing pixel observations. Since locations outside the field of view will have valid observations in at most two cardinal directions, pixels with valid values in at least three cardinal directions were filled with nearest-neighbor interpolation; qualitatively, this better resembled level-1 data and did not assign improbable ”cold” pixels islands in the middle of large lava lakes/flows.

To narrow the focus of the image forecasting task, a $96\times 96$ pixel subset of the scene was used in the data set; with a 90-meter spatial resolution, this translates to an $8.64\times 8.64$ kilometer subset. This precise size is somewhat arbitrary, but, qualitatively, it was large enough to observe significant portions of large lava flows while not being too large to neglect small lava lakes. The subset was centered on the latitude and longitude of the volcano summit. Converting the scenes’ coordinate reference systems to latitude and longitude for summit locating and subsetting introduced geometric distortions that were visually apparent for Mount Erebus (located in Antarctica); however, these relics were not problematic for thermal interpretation.

In volcanology studies, temperature measurements are often examined relative to the environmental background (i.e., the excessive temperature of the volcanic feature relative to the surrounding environment); as such, the AST_08 scenes were processed this way (see Figure 3 for steps up to this point). Volcano studies often manually identify a cloud-free, representative part of the environment to determine the environmental background temperature. To automate this process in a robust way, the following process was used to compute the average environmental background temperature: first, a $10\times 10$ pixel subset was extracted from each corner of the reduced image; these corners were then considered appropriate if less than 10% of their pixel values were outside the sensor’s field of view (i.e., missing values). The mean was then calculated for each subset, and the mean of those means was used as the environmental background temperature. Mount Erebus’s location along the orbital path led to unique viewing geometries where this method would fail to identify any valid corners; in those cases, the $10\times 10$ pixel subsets were rotated around the perimeter of the image to identify suitable areas (see Figure 4). Once the environmental background temperature was computed, it was subtracted from every pixel value. Notably, in addition to being consistent with volcanology studies, this removes the task of forecasting the seasonality of the background environment.

After the above, there remained missing values from the $96\times 96$ scene subsets that represented areas outside the sensor’s field of view. These pixels were then filled with the values from the previous scene and the date of those previous scenes was noted. When this occurred for the first scene acquired for a volcano (i.e., no previous scenes), the missing values were filled with nearest-neighbor interpolation (Figure 5). Lastly, prior to training the data was scaled between 0 and 1 (i.e., $x_{new}=(x_{original}-x_{min})/(x_{max}-x_{min})$) to facilitate training convergence.

The forecasting task here uses the previous $n$ scenes to forecast the next scene. For most models, this is done by treating each pixel as an independent observation, but for convolution-based neural networks, this is done at the image-level. The exact choice of $n$ is not straightforward, but realistically, it must be large enough to make use of past information and small enough to not magnify computational and storage burdens. A popular method in volcanology studies is to examine the maximum excess temperature through time, so these time series were extracted from the L1T data (prior to AST_08 retrieval) and their autocorrelation structure was examined for each volcano. The average of the largest significant autocorrelation lag was $6.43$, so $n$ was set to $6$.

In the case of volcano-specific training sets (see Results and Discussion, section C), the largest significant autocorrelation lag for each volcano was used with a minimum accepted value of $3$ and a maximum accepted value of $10$ previous scenes.

The data set was separated with a $70/15/15$ split chronologically by each volcano into a training, validation, and test set. The training set was used to fit model parameters, the validation set was used to compare different models as well as compare the same model fit with different regularization strengths, and the test set was aside for evaluation of the final, best model. Every model was fit with $5$ different regularization strengths ranging from $0.0001$ to $1$ by a factor of $10$; higher values of this parameter help reduce overfitting to the training set by penalizing large model parameters.

Choosing the exact number of hidden layers as well as the number of units/neurons in each layer is not straightforward, with the only general rules being that more parameters increase performance, computational burden, and the need to combat overfitting. The exact hidden layer choices used in this project were taken from the published best results of the existing neural networks used here [51, 55, 56]. All models were trained for 100 epochs with the Adam optimizer [59] to minimize mean squared error. Optimization occurred in varying mini-batch sizes that were determined to maximize GPU utilization.

All models were implemented in the Python programming language via the PyTorch machine learning framework. The hardware used for this project was an HPC cluster with 10 NVIDIA Tesla V100 GPUs with 32 gigabytes of RAM each; all neural network computations used this hardware in parallel.

Table 1 presents the validation set RMSE along with important training details for all the models. Among the naive methods, using the most recent scene as the forecast for the next scene led to an RMSE of $5.291^{\circ}$C while, alternatively, assuming a homogenous temperature landscape (“All Zeros”) led to an RMSE of $4.738^{\circ}$C. Representing the only learned statistical model, the AR(6) returned an RMSE of $4.322^{\circ}$C. Among the non-convolutional neural networks, the standard LSTM models returned RMSEs ranging from $4.217$ to $4.265^{\circ}$C, the Time-LSTM models returned RMSEs ranging from $4.694$ to $4.870^{\circ}$C, and the Time-Aware LSTM models returned RMSEs ranging from $4.660$ to $4.726^{\circ}$C. The latter two converged very closely to the All Zeros model, returning low magnitude predictions in the best results. Convolution-based neural networks returned RMSEs between $4.164$ and $4.428^{\circ}$C with less variance than the non-convolutional neural networks. The lowest RMSE was achieved by the proposed ConvLSTM + Time-LSTM + U-Net architecture. It was approximately $10$ times faster than the lowest RMSE LSTM and required approximately $40$ times fewer parameters.

For the best model of each type (i.e., statistical, non-convolutional, and convolutional), the models were separately trained for $1000$ epochs to determine if further training resulted in better performance, but, instead, it was found that validation set RMSE increased for all models with this extended training. For the AR(6), validation set RMSE increased from $4.322$ to $4.357^{\circ}$C; for the LSTM, $4.217$ to $4.293^{\circ}$C; and, for the ConvLSTM + Time-LSTM + U-Net, $4.164$ to $4.182^{\circ}$C. The models originally trained for $100$ epochs are thus used for subsequent discussion, and the three mentioned in this paragraph will be referred to as the “selected models”.

Predictions for the lowest RMSE model are presented in Figure 6 and Figure 7. Both of these demonstrations compare observations with predictions during large effusive events. In Figure 6 the predictions are well aligned with the observations. It can be noted that the predicted sequence failed to predict new spatial features until after they emerged in observations. Interestingly, as the lava flow continues, the model emphasized particular hotspot areas that were present in prior scenes while broadening the width of more linear features. While Figure 7 is also a large lava flow, its data are much less ideal—note the persistent presence of partial cloud cover as well as missing edge values which required filling with previous scenes. Despite these obstructing conditions, the model continued to predict the presence of a large lava flow in the correct general direction and shape.

The distributions of predicted pixel values for the selected models as well as a poor-performing model (Time-LSTM) were examined and are presented in Figure 8. Barring the poor-performing model, it was found that the models with lower RMSE tended to predict lower-valued extremes than higher RMSE models and all models predicted lower-valued extremes than the true observations.

All in all, the explicit modeling/optimization task was to predict image sequences that are irregularly spaced in time with minimal RMSE relative to the true observations. Each model had different characteristics which influenced its ability to accomplish this task. For example, the AR(6) had very few parameters while some of the neural networks had millions of parameters; the Time-LSTM considered the time between scenes while the ConvLSTM considered a pixel’s spatial neighbors. Ultimately, the proposed model, specifically the ConvLSTM + Time-LSTM + U-net, considered more characteristics of the data (e.g., spatial, irregular in time, and high resolution) and accomplished this explicit goal the best. It is worthwhile to note that while all models have low error relative to the full range of the data that these error differences are spread across $2211840$ pixels—$240$ images in the withheld validation set, each of which contains $9216$ pixels.

With infrared-based imagery, there are many derived time series of interest in volcanic monitoring. Specifically, with temperature data, radiance-based measurements cannot be examined, so here we focus on (1) maximum temperature above the background, (2) the number of hotspots, and (3) distance of furthest hotspot from the summit. All derived time series are concerned with anomalous or extreme temperatures, so it was expected that the models’ inability to predict the true scale of extreme temperatures would inhibit derived time series. To alleviate this, predicted time series were derived both before and after histogram matching the validation set predictions to the training set observations. A demonstration of this can be seen in Figure 9.

Table 2 contains the error of the selected models as well as the poor-performing Time-LSTM on the derived time series. See Figure 10 for a demonstration of these time series using the validation set, histogram matching, and the ConvLSTM + Time-LSTM + U-Net model.

Before histogram matching, the AR(6) achieves the lowest RMSE on all three derived time series. This is attributable to the fact that it predicts the largest maximum values of the available models. After histogram matching, which accounts for the apparent tendency of neural networks to produce low-scale extremes, the AR(6) model is no longer the best at any derived time series. Likewise, the proposed model (ConvLSTM + Time-LSTM + U-net), which achieved the lowest RMSE on predicting the image sequences was unable to reproduce any derived time series with the lowest RMSE. Instead, after histogram matching, the traditional LSTM achieved the lowest RMSE recreating the maximum temperature above the background ($32.120^{\circ}$C), while the Time-LSTM achieved the lowest RMSE on the number of hotspots ($17.468$ pixels) and the maximum distance of a hotspot from the summit ($825.800$ meters).

This discrepancy is interesting because the Time-LSTM was included in further comparisons since its best optimization resulted in the worst validation set image sequences (RMSE = $4.694^{\circ}$C). Further examination revealed that this is because it tended to converge to near-zero image predictions (minimum = $-0.907^{\circ}$C, maximum = $0.958^{\circ}$C). Together, this suggests that despite struggling immensely with scale, the Time-LSTM learned superior pattern recognition within a very narrow comfort zone for this task that was indirectly optimized for.

Scientifically, this suggests an important implication—when forecasting multipurpose data with many potential uses, directly minimizing error on the data does not necessarily minimize error on its potential uses or derivations. Moving forward, it would be best to explicitly account for these other metrics directly in the loss function (i.e., optimization problem). Compared to highly engineered, domain-based algorithms, deep learning models are relatively simple and rely on heavy parameterization along with large amounts of data to improve performance. With immense interest in the value of deep learning for scientific tasks, it will likely be a very promising pursuit to engineer loss functions with the multiple scientific priorities and constraints in mind—an example of this is can be seen when forecasting water temperature constrained by physical laws [60].

Broadly, RMSE decreased when the training set was limited to individual volcanoes (Table 3). In general, volcanoes with fewer data produced worse results, but performance did not increase exclusively with training set quantity. For example, training exclusively with Erta Ale or Kīlauea data led to lower RMSE than training exclusively with Erebus data ($4.356^{\circ}$C, n = 178 and $4.390^{\circ}$C, n = 150 versus $4.684^{\circ}$C, n = 336), and training exclusively with Nyiragongo data led to higher RMSE than training exclusively with Masaya data ($5.365^{\circ}$C, n = 74 versus $5.343^{\circ}$C, n = 44).

It is commonly noted that every volcano is different and, as such, patterns learned from one volcano may not generalize to another. It could have been possible that a volcano’s validation set performance was minimized by training only with data from that volcano, but that was not found to be the case; in general, training with a specific volcano did not reduce RMSE on that volcano itself. The only exception to this was Erta Ale (all volcanoes RMSE = $5.875^{\circ}$C, Erta Ale only RMSE = $5.534^{\circ}$C). The reasons for this are not known; it may be that patterns from the other volcanoes do not apply well to Erta Ale, and/or the fact that, unlike most volcanoes in this data set, Erta Ale has relatively flat, cloud-free terrain.

In $8/9$ cases, a training set consisting of one volcano did not benefit validation set performance on that volcano itself, but in $4/9$ cases, this did reduce validation set RMSE on at least one other volcano. This was the case when the training set consisted of only Erebus, Erta Ale, Kīlauea, and Pu’u ’Ō’ō resulting in validation set RMSE reductions in $2/8$, $3/8$, $5/8$, and $1/8$ other volcanoes, respectively. Among these, the model fit using only Kīlauea data improved validation set performance on at least the same volcanoes that the others did (i.e., training sets of Erebus, Erta Ale, and Pu’u ’Ō’ō). Kīlauea is one of the most well-studied volcanoes in the world and, fortunately, this appears to suggest that thermal patterns learned from it may generalize better than other well-known basaltic volcanoes.

A forefront weakness in this forecasting project is that high-resolution thermal imagery is the only data considered. Even within thermal data sources, there are additional data sources that could be integrated—most obviously MODIS. Additionally, regional weather satellites such as GOES could prove valuable and would still provide nighttime imagery. Beyond nighttime imagery, there are also the options of Landsat and Sentinel-2. With additional expertise it would also be beneficial to include other forms of remote sensing, such as InSAR for ground deformation and IR/UV for gas emissions, and, when available, performance would likely benefit from the inclusion of ground monitoring efforts. Ultimately, volcanoes display a variety of behavior that can be measured on the ground and from space; here, we focused on thermal satellite imagery, one of the most globally available and widely adopted.

Moving forward, the task of integrating all of this data will likely not be straightforward. Experts in the field have speculated that ”few developments promise to revolutionize eruption forecasting — along with so many other branches of science — as does machine learning.” [1, p. 23]. We argue here that important developments moving forward will likely require the development of machine/deep learning methods with scientific goals and data characteristics in mind. Geoscience fields, for example, are generally data-sparse but rich with theory-based models [61], and it has been found that integrating scientific theory with machine learning helps alleviate the burden of low data quantity [60]. Additionally, scientific data may not have the same characteristics as the data that machine learning methods were developed for. For example, Earth data are often irregular in time and space, and imagery may be very different (e.g., gravity-based data, interferograms, hyperspectral) from the typical optical situations that computer vision is concerned with. In summary, there are many broad directions and considerations for machine learning in science, and scientific discovery and application will likely benefit from models that are creatively specified to not only the data but also the governing scientific goals and principles.