Volcanic smoke reduction in visible and thermal infrared imagery

Abstract

A simple and rapid technique to effectively remove volcanic smoke (vapor-dominated white plume) from visible or thermal infrared images of active craters is proposed, in which multiple photographs taken in a fixed frame are post-processed by a script-based program on a computer. This method allows us to obtain a detailed view of the surface conditions on fuming craters that is often impossible to observe by a single snapshot. It is also effective in improving an accuracy of surface temperature anomalies, and hence, heat flux from hot lakes or geothermal areas.

Introduction

Photography is one of the fundamental methods of visual observation in volcanology. Currently, digital cameras are in common use and are occasionally affixed to an unmanned aerial vehicle (UAV), allowing visible and thermal infrared (TIR) imagery to be obtained for accurate surveillance of crater morphology, thermal anomalies, photogrammetry, etc. Obviously, the greatest benefit of digital photographing is affinity to automatic recording and computer-based processing. For instance, Terada et al. (2008) monitored the water level of the acidic crater lake of Aso volcano for 15 months by automatic interval recording of photographs from a fixed ground site. Rouwet (2011) collected and analyzed the photographs that were manually taken from the same spot at the crater rim of El Chichón volcano to discuss the morphological changes of the crater lake during the period 2001–2007. Harvey et al. (2016) used a drone with a TIR camera to obtain a high-resolution georeferenced thermal imagery of the Waikite geothermal area, New Zealand. Moussallam et al. (2016) constructed a 3D model of the crater morphology by means of photogrammetry technique using the video footage from a UAV that flew above the lava lake at Villarrica volcano. These are successful results under favorable conditions. However, in our experience, it easily becomes difficult to obtain a clear view of an active crater when volcanic smoke conceals the site. In many cases, the crater bottom is covered with white smoke, partly hindering the observer from obtaining an overall view of the surface details. Thus, a single snapshot taken from the crater rim does not always provide us much useful information. Nonetheless, as long as we have some windows through the smoke as it shifts due to variable winds, it seems possible to obtain a composite smokeless image by gathering the visible parts into a single picture as a mosaic. One of the conventional ways of doing this is hand sketching. Today, the same procedure can be easily accomplished by digital compositing, in which a computer performs the posterior image processing. In this case, field observations involve taking multiple snapshots or a video image from a fixed position for a certain time period. A benefit of using digital imagery instead of hand sketching is that we can easily obtain a high-resolution objective record regardless of the sketching skill of the observer.

Herein, we propose a simple and rapid method to effectively remove the volcanic smoke component by compositing multiple photographs that are taken from a fixed position. This image processing technique enables us to recognize the details of a concealed crater wall and/or bottom (e.g., the number and size of fumaroles, existence of sublimed sulfur, thermal ponds). Moreover, application of this technique to TIR images is also useful for removing the effect of relatively cold and weak plumes that are generally semi-opaque to infrared, which improves the accuracy of estimates on heat discharge rate from ground surface or hot crater lakes. We show the algorithm and sample code of the image processing as well as the examples of field application at some volcanoes.

The algorithm of smoke reduction

Smoke without solid particles such as volcanic ash generally looks white in visible imagery, and thus in many cases, brighter than the background wall rocks, vegetation, soil, or sky. Looking at a specific position on the time-sequenced photographs, the brightness may change with time. It is expected that the brightness will be relatively high at the moment when the position is occupied by fume. In other words, the lower the brightness is, the less significant is the influence of the smoke in the line of sight. Based on this concept, image processing is carried out using the following procedure.

As the first step, suppose we have N sheets of photographs, which are not necessarily taken at a constant time interval. We name the first and second images as X and Y, respectively. We compare the brightness of the same pixels on X and Y and choose the one with lower brightness. This procedure is performed throughout the images of X and Y. A composite image Z is constructed using the chosen (darker) pixels. Second, the composite image Z is renamed as the new X, and the third photograph is now called Y. Then, the comparison procedure repeats for the new X and Y, replacing the composite image Z. By repeating the procedure N − 1 times, we obtain a final composite image that is composed of the lowest brightness for each pixel among N photographs.

In our processing script, we define the brightness of a pixel A_ij at the position (i, j) as,

$$A_{ij} = \frac{{R_{ij} }}{3} + \frac{{G_{ij} }}{3} + \frac{{B_{ij} }}{3}$$

(1)

where R, G, and B represent the intensity of red, green, and blue light, respectively. Note that replacing the weighting of RGB with that of luminance, namely (0.299, 0.587, 0.114) instead of (1/3, 1/3, 1/3), makes no significant difference. Comparison and selection of the scene of lower brightness can be realized by performing the following operation,

$$Z_{ij} = \left( {A_{Xij} < A_{Yij} } \right) \times X_{ij} + \left( {A_{Xij} \ge A_{Yij} } \right) \times Y_{ij}$$

(2)

where inequalities in the parentheses represent the comparison operation and their results are supposed to be given as a Boolean value (1 for true, 0 for false). X_ij, Y_ij, and Z_ij are the RGB data for each pixel of the images X, Y, and Z, respectively. The above procedure is schematically summarized in Fig. 1. The core of the compositing procedure is only composed of these two operations and can be written in a very simple and short script, which is available in the additional files of this article. In short, we gather the minimum brightness for each pixel and composite into one picture.

Fig. 1

Schematic flow of the compositing image processing of this study

Before we show the applications to volcanoes, we illustrate how this composite processing works using the cloud images obtained from Himawari-8, a meteorological satellite, as a demonstrative example. In the cloud images from the satellite, brightness is higher (whiter) in areas with a significant cloud cover, from which stronger IR light is radiated, whereas bare land and sea are colored in green and blue (namely, darker than clouds), respectively. Clouds are advected by wind, and their shapes are deformed with time. Such features of the cloud image are also characteristic of the photographs taken in volcanic fields. The cloud image is provided by the Japan Meteorological Agency every 30 min on their Web site, and it is suitable for verifying the performance of our image processing method. We show the cloud images at 15:00 (UT) August 9, 2015, and 15:00 (UT) on August 11, 2015, respectively, in Fig. 2a, b. The result of compositing these two images by our algorithm is presented in Fig. 2c. In the composite image, it can be seen that clouds are partially removed by selecting the pixels with a relatively low brightness. Figure 2d shows the result of applying this process to 95 time-sequenced images, taken at a time interval of 30 min. Clouds are now effectively removed throughout the image except for thin clouds remaining in the inland region of China, part of the Pacific Ocean in the southeast of the Japanese archipelago, and near the Kuril Islands located in the northeast of Hokkaido.

Fig. 2

Image data are downloaded from the Web site of Japan Meteorological Agency (http://www.jma.go.jp/jp/gms/smallc.html)

Satellite IR imagery of Himawari-8. a IR image at 15:00 (UT) August 9, 2015. b IR image at 15:00 (UT) August 11, 2015. c Composited result of the two images a and b. d Composited result of 95 images from August 9 through August 12.

Although some commercial image processing software is equipped with similar functions, including compositing as one of the variety of retouching effects, they are normally operated on a one-by-one basis and hence require considerable time and effort to deal with multiple images. In contrast, our method is script-based and is suited to dealing with multiple images in an automatic way. In our current environment (MATLAB^® R2016a, on Apple Macbook Air 1.8 GHz Intel Core i5) processing time is approximately 180 s for 50 images of 6016 × 4000 pixels and 15 s for 60 images of 1600 × 1200 pixels.

Applications to observation of volcanic craters

62-II crater, Mt. Tokachidake

We first show an example of photographs taken from the rim of 62-II crater, one of the active craters of Mt. Tokachidake volcano, central part of Hokkaido, northern Japan on July 8, 2013. The photographing interval is not constant, but occurred at roughly 5–30 s. The quasi-circular crater measures roughly 150 m in diameter. Figure 3 compares an original snapshot from the north wall and the result of smoke reduction using 60 shots. In the original image, the pit craters are covered with white smoke and it is not possible to recognize the number and size of the pits in detail. In contrast, in the composite image, the smoke is significantly removed, and it is clearly seen that sublimed sulfur adheres to the fumarolic vents. More specifically, in the central part of the screen, we can also recognize the number and size of pit craters that look like dark spots, probably indicating that the sublimed sulfur turns dark due to high temperatures at the root of the fumaroles on the crater wall.

Fig. 3

An example of application to 62-II crater, Mt. Tokachidake, Japan. a A snapshot before processing. b Composited image of 60 snapshots. The white bar in the right panel indicates an approximate scale of 10 m

Naka-machineshiri crater, Mt. Meakandake

As an another example we show the pictures taken from Naka-machineshiri crater of Mt. Meakandake volcano, eastern Hokkaido. As the crater hosts numerous fumaroles and small hot ponds on its 250 m wide bottom, it is difficult to obtain a clear view of the whole crater from any single snapshot. Figure 4a is one of the original images prior to processing that was taken on June 25, 2013. Since white smoke widely covers the central part of the image, it is impossible to observe the crater’s bottom in detail. However, composite processing using 50 images in about 10 min effectively removes the smoke as shown in Fig. 4b. Enlarging the central part of the processed image (white rectangle in Fig. 4b), we recognize separate vents (No. 1–7 in Fig. 5) of several to 10 m wide at the roots of the plumes, some of which are hot ponds that are fountaining (Fig. 5, No. 4). In addition, it can be seen in the right side of the close-up image that hot spring water is running down as a creek toward the valley at the back center of the image. Downsized sample photographs and a MATLAB^® script are available in Additional file 1.

Fig. 4

An example of application to Naka-machineshiri crater, Mt. Meakandake, Japan. a A snapshot before processing. b Composited image of 50 shots. Black dots indicated with arrows are birds. The white bar in the right panel indicates an approximate scale of 20 m

Fig. 5

A close-up view of the crater bottom (the white rectangle in Fig. 4b). Numbers 1–7 indicate pits or hot ponds. Fountaining is seen at the hot pond No. 4. The white bar in the panel indicates an approximate scale of 20 m

Applicability and precautions

Our method of smoke reduction processing can be used for various applications in volcanology. One of the most basic usages is visual observation and monitoring of superficial low-intensity volcanic activity. As illustrated in the previous section, detailed description of the number and size of fumaroles and hot ponds may be useful in estimating the heat flux. Morphology inside a crater may also be valuable information. Our method, when combined with the recently developed technique of digital photogrammetry (e.g., Lewis et al. 2015; Thiele et al. 2017), will enable us to obtain a digital elevation model of a fumed crater in the absence of interfering smoke. Further, applying this method before and after an eruption will be effective in estimating the amount of ejected deposits and changes in topography by the eruption.

However, the following points are notable. Flying objects, such as birds or aircrafts, can remain as dark dots in the final composite image, as they are visually darker than the background sky (see Fig. 4b). Another concern is that the overall tone of the final composite image is prone to be darker than individual raw images, as this method always prefers a scene with relatively lower brightness. Therefore, care must be taken when using the images captured under such conditions that intensity of sunlight changes significantly during the period of observation, or that patches of clouds in the sky are shifting in front of the sun. In addition, this method is unsuitable for removing smoke photographed with backlighting or a plume that contains dense ash, since it is based on the premise that the smoke or plume will look brighter than the background.

Next, we investigate how many images are necessary to effectively remove the smoke. To this end, we examine the previous examples of Mt. Tokachidake and Mt. Meakandake more closely. The initial image before composite processing is regarded as the state in which the smoke is not removed at all, whereas the final composite image (60 shots for Tokachidake and 50 shots for Meakandake) is regarded as the state in which the smoke has been fully removed. Here we define the degree of smoke reduction (DSR) as an index indicating the extent of smoke removal after the i-th compositing step relative to the initial image,

$${\text{DSR}} = N_{\text{r}} \left( i \right) \, /N_{\text{r}} \left( {\text{imax}} \right)$$

(3)

where N_r(i) represents the number of pixels in the i-th composited image that experienced the brightness reduction relative to the initial one by 5% or more. It should be noted that N_r depends on i, the number of compositing, and that N_r(imax) is given for the final composite image. Figure 6 depicts the change in DSR as a function of i. By definition, DSR in the initial state is equal to zero, while that in the final state is 1.0. As can be seen in Fig. 6, DSRs calculated in the order of time (thick lines with dots) show some stepwise increase in the course of compositing, reflecting the sudden clearance of view that had previously been concealed with smoke at the time of change in the wind direction. In either the case of Tokachidake and Meakandake, DSR exceeds 0.9 at around i = 30; thereafter, DSR is only improved in a gradual manner.

Fig. 6

DSR (see text for details) with respect to the number of pictures composited. Blue and red curves show the examples from Mt. Tokachidake (Fig. 3) and those from Mt. Meakandake (Fig. 4), respectively. The thicker lines with symbols are processed in the order of time, while the thinner lines are the cases in which pictures are randomly resorted and composited, and averaged over 10 times

We also investigate the cases in which the images are composited in a random order with respect to time. Here, composite processing with the randomly sorted images is repeated ten times and then averaged. The result is shown as two thin lines in Fig. 6, in which previous sharp steps are now smoothed out. Further, it can be seen that DSR reaches 0.8 around i = 30 and 0.9 around i = 40. Although the number of sheets necessary to obtain a satisfactory effect generally depends on the amount of smoke coverage, wind conditions, interval of photographing, and so on, we recommend i = 40 as a guide, for which the interval of observation is several seconds to some tens of seconds. Reduction of shooting interval does not always result in effective smoke reduction when the smoke is shifting slowly relative to its size scale.

Applications to TIR imagery

Crater lake, Aso volcano

Spampinato et al. (2011) provided a critical review of the application of TIR cameras to volcano surveillance and monitoring, in which they discussed several potential error sources that can influence derived temperatures, such as emissivity of targets, atmospheric attenuation, viewing angle, instrumental characteristics. These factors are sometimes significant and cannot be neglected. Aside from these errors, white plumes in front of the targets are the primary nuisance for TIR imagery in crater areas or steaming grounds, since they are semi-opaque to infrared. Applying the compositing algorithm of this study to TIR imagery greatly reduces the problem. Here we provide an example using the hot crater lake of Mt. Aso volcano, southwest Japan.

The first crater of Aso volcano hosts a hot acidic crater lake, which is locally called Yudamari in non-eruptive periods (Terada et al. 2008). Temperature of the crater lake provides important information in evaluating the thermal discharge rate from the lake surface (Terada et al. 2012; Terada and Hashimoto 2017). Figure 7a shows a snapshot of the surface temperature of Yudamari lake that was taken by an TIR camera on February 2, 2009. Composite processing of 120 scenes by means of the method in this study yields Fig. 7b. A visible image with an approximately equal frame size is shown in Fig. 7c for reference. The lake had a diameter of approximately 200 m at the time when it was photographed. The observation was performed from the western crater rim of the first crater in about 30 min. Distance from the lake center to the observation point was roughly 300 m, and the relative height from the lake level was about 100 m. We note that in this case we did not evaluate the color data of the images, but temperature values were calculated in a comparison operation as shown in Eq. (4) instead of Eq. (2),

$$Z_{ij} = \left( {T_{Xij} < T_{Yij} } \right) \times X_{ij} + \left( {T_{Xij} \ge T_{Yij} } \right) \times Y_{ij}$$

(4)

where T_Xij and T_Yij indicate the apparent temperature of the pixel position (i, j) of the images X and Y, respectively. Histograms of the apparent temperature for the TIR images are shown in Fig. 8, in which we recognize a significant rightward shift and convergence of high-temperature component, while the background component of less than 20 °C stays almost unchanged. Conversion from surface temperature to heat discharge using the method of Terada et al. (2008), based on Ryan et al. (1974), results in an increase of approximately 20%, from ~ 18 to ~ 22 MW.

Fig. 7

Crater lake of Aso volcano. a An IR thermograph before processing. b A composited TIR image using 120 scenes. Assumed emissivity is 1.0 (instrumental default). No atmospheric corrections applied. c A visible image of the crater lake. The white bar in the panel indicates an approximate scale of 20 m

Fig. 8

Histograms of the apparent temperature in the two images a and b of Fig. 7. The histogram of a is superimposed as transparent bars on panel b for reference

As is seen from the visible imagery (Fig. 7c), evaporated water from the lake floats above the lake surface as white steam and is gently convected by the wind inside the crater, forming some distinct filaments. It can be seen from the TIR image prior to the composite processing (Fig. 7a) that the apparent temperature on such vapor filaments is considerably lower than other areas. High-temperature domains, roughly tens of meters across, are divided by these filaments. Low-temperature filaments have been removed in the composited image, rendering the surface temperature distribution quite smooth (Fig. 7b). Temperature zoning is also recognized, in which the apparent temperature gradually decreases from bottom (~50 °C) to top (~43 °C) of the image. Such zoning of the apparent temperature may be attributed to the atmospheric attenuation (Spampinato et al. 2011) by distance between the observation point and the object, which is mainly due to unsaturated (transparent) vapor in the air along each line of sight. Namely, the greater the distance is, the more intense the attenuation is. Accordingly, it is likely that the apparent temperature of the lake surface in the right bottom side of the lake on the image (47–48 °C) is closest to the actual water temperature. Given that the IR thermometry gives systematically lower values than direct measurements from a buoy by 8.3 °C, which is empirically confirmed by Terada et al. (2012), and assuming the skin effect (slight decrease in water temperature in a boundary layer at the very surface of the lake; Oppenheimer 1997) as 2 °C (Terada et al. 2008), the real water temperature on that day is estimated to be 53–54 °C. Besides the zonal distribution, we recognize several hot spots in the central part of the lake. Undulation due to the steam filaments mentioned above does not account for these hot spots, since it is confirmed that these spots are stable regardless of additional compositing. We consider that these hot spots correspond to the upwelling of high-temperature water or steam that has been envisaged based on the heat and water balance of Yudamari by Terada et al. (2012). The image processing of the present study can contribute to the detection of such spot-like temperature anomalies also for other crater lakes.

Steaming ground, Mt. Tokachidake

There are numerous studies on the heat radiation from geothermal anomalies and steaming grounds using IR thermometry (e.g., Sekioka and Yuhara 1974; Stevenson and Varley 2008; Harris et al. 2009; Spampinato et al. 2011). Accurate mapping of thermal anomalies is of the fundamental importance in heat radiation estimation, where apparent reduction of temperature due to water vapor is a major source of estimation error. We show here in Fig. 9 an example from a steaming ground at Mt. Tokachidake. Thermal infrared images were acquired on the steaming ground situated on the southern side of 62-II crater. Our viewpoint was from 1 km south, across a valley, where we had an overall picture of the steaming ground within the viewing angle of the TIR camera. The observation was conducted in the early morning on July 30, 2015, in order to avoid the effect of solar radiation. We took 26 shots in approximately 35 min. By compositing the TIR images in the same way as described in the previous section, we obtained a composited thermal image in which the thermal anomalies were enhanced. The actual area for a single pixel corresponds approximately to 0.67 m². A histogram of the apparent temperature is presented in Fig. 10. Percentage of pixels over 25 °C in the images before and after the composite are 4.5 and 6.2%, respectively. The heat flux is estimated based on the method of Sekioka and Yuhara (1974) as ~ 7.9 and ~ 9.2 MW for the raw and composited images, respectively. We should note that the absolute values of the heat flux may be somewhat underestimated as the distance from the observation point to the target area is large, and thus, suffers from atmosphreric attenuation (e.g., Spampinato et al. 2011). Nonetheless, we consider the relative improvement of heat flux values by means of composite processing to be significant and meaningful for monitoring purposes, as long as the atmospheric conditions can be ragarded as constant.

Fig. 9

A steaming ground on the southern side of 62-II crater, Mt. Tokachidake. a An IR thermograph before processing. b A composited TIR image using 26 sheets. Assumed emissivity is 1.0 (instrumental default). No atmospheric corrections applied. The white bar in the right panel indicates an approximate scale of 100 m. c A photograph whose frame roughly corresponds to the TIR images

Fig. 10

Histograms of the apparent temperature in the two images of Fig. 9. The histogram of a is superimposed as transparent bars on panel b for reference

Conclusions

We proposed a simple and convenient method of removing volcanic smoke from multiple photographs taken in a fixed frame. Applying this method to visible images enabled us to observe the details inside fuming craters. We confirmed from some examples that compositing of approximately 40 shots results in satisfactory smoke reduction, although it intrinsically depends on how densely smoke conceals the site. Application of this method to TIR imagery greatly reduced the effects of IR absorption due to water vapor or mist in the smoke, which helped us to determine the surface temperature of crater lakes or geothermal areas more accurately and to improve the accuracy of heat flux estimation. We note that our method is also applicable to emerging observational platforms such as UAV-based imagery and photogrammetry as long as the vehicle can hover at a suitable position in the air for a few minutes.