1. Introduction
My mind tends to bring this visual to the forefront when I encounter the word Volcano-
Those who may not be familiar, Sabu is a giant alien who arrived on Earth from his home planet Jupiter and resides with the diminutive yet clever old man, Chacha Chaudhary - also the name of the very popular comic book series in India. Together, they form a perfect pair of brawn and wit who help the society by getting rid of notorious elements. While being proficient at brushing aside crooks, Sabu does tend to struggle against a scheming or occasionally, an even stronger adversary. When the lives of Chacha and his family are threatened, Sabu enters into a Hulk-esque rage which increases his strength manifold and is a forewarning about the destruction that will ensue. Pran, the creator of the comic series, signals this to the readers using the erupting volcano depiction as seen in Figure 1.
A big thank you to RUS Copernicus & EO College - I have utilized their learning resources, as I have done on several occasions in the past, to form an understanding about this topic and to create the demonstration that follows in this post.
SECTION HYPERLINKS
2. Setting the Context
Humans have feared and revered volcanoes through time. While India is home to just a couple of active ones, far away from the mainland, we do have a striking name for volcanoes - Jwalamukhi - which can be literally translated from Hindi as the mouth that spurts fire.
Much like a human mouth but at a longer timescale, a volcano breathes in and out i.e. its size, shape and volume expands and contracts due to the changes in pressure exerted by the hot magma within, which the surrounding rocks accommodate by deforming. The trend of Deformation is what volcanologists observe keenly as it is a leading indicator of eruption potential - the more restless the magmatic system becomes, it causes an increase in both, the frequency as well as the intensity of 'uplift' (ground level rising) and 'subsidence' (ground level sinking). Thus, volcanic unrest causes ground to displace and deform, which often is a precursor to an eruption event.
Researchers use geodetic monitoring techniques such as Leveling, GNSS & InSAR to identify and measure Deformation - it is not discernible to the naked eye.
In this post, I will demonstrate the method of generating a Deformation Map and deriving the quantitative rate of Displacement using the latter technique - Interferometry (In) applied on Synthetic Aperture Radar (SAR) Satellite Imagery.
As you continue reading this post, you'll find that I've used InSAR to monitor Deformation at three recent Eruption sites-
Monitoring Deformation at a Subsurface level i.e. close to the source of volcanic activity - the Magma Chamber which lies 1-10 kilometers below the Earth's surface - would definitely be more meaningful than monitoring Deformation at/above Ground Level near the eruption outlet, just that the current technology is not able to penetrate the kilometers of rocks that lie between the surface and the underground chamber😓.
The difficulty is compounded when there is an additional layer of water in between. Our inability to monitor Deformation in submarine conditions is worrisome as this is where most of our planet's productive volcanic systems lie - we had no inkling that Hunga Tonga–Hunga Haʻapai was to erupt in a violent expression of fury on 15th January 2022.
3. The Processing Chain and Video Walkthrough
As indicated earlier, I will apply Interferometry on Sentinel-1 SAR Imagery using SNAP tool to map Ground/Surface Deformation and to derive the rates of Displacement (Uplift and Subsidence).
The Processing Chain involved in doing so is extensive, here's a graphical representation-
I have attached a step-by-step Video walkthrough of this workflow below, the area of study being Cumbre Vieja Volcano in Canary Islands, Africa which erupted in 2021. Those who prefer to watch content over reading it will stand to benefit from this visual guide. However, I'll recommend that you read this post and choose to view this upon completion.
Video 1: Deriving Deformation & Displacement Maps for Volcanic Eruption (Cumbre Vieja 2021) using Synthetic Aperture Radar Satellite Imagery
4. Step-by-Step Guide on InSAR Interferometry to map Deformation
Let me begin to explain the InSAR Processing Chain as depicted in Figure 3. Download the image to your PC for easy reference (left-click to expand the image and then right-click it to save it).
PS: In the visual that accompanies each section note below, I often alternate between the three eruption events I have chosen to use in this post - the image description specifies the study area
a) Read & Read(2) - Two sets of SAR imagery are used as the input datasets in the workflow - one acquired prior to the eruption and the other after the eruption event. Essentially, I am asking SNAP to read these input datasets which I will subsequently compare / perform change detection, in order to detect Deformation and measure the rate of Displacement.
b) TOPSAR-Split - Satellite Imagery datasets tend to be large - each Sentinel-1 SAR dataset that I am using is around 1 GB in size and has a coverage (swath width) of 250 kilometers and spatial resolution is 5 x 20 m per pixel. The TOP in TOPSAR refers to Terrain Observation with Progressive scans - a method of illuminating targets which Sentinel-1 satellite uses, and which is particularly useful for Interferometry and Time Series analysis.
In order to reduce computational load and save time during subsequent processing, I would like to restrict the spatial extent of the Imagery dataset to just the the region of my interest - the volcanic eruption site, as well as restrict the number of spectral bands containing reflectance information to just the one which I need for this demonstration - VV polarization. The TOPSAR-Split operator helps me to accomplish these objectives of spatial and spectral subsetting.
c) Apply Orbit Files - A Satellite captures data as it navigates in its Orbit. The Orbit Files within the SAR Imagery dataset package contain the satellite's velocity and position metadata at the time of Imagery acquisition - this radar geometry information helps in the precise assignment of energy reflectance values that the satellite captures to the individual pixels in the Imagery dataset. The reason why it is necessary to perform this step is because refined Orbit Files take up to a few weeks to be prepared by the satellite operator and thus, the original ones bundled in the Imagery Package may not be entirely accurate, particularly if the Imagery has been acquired very recently, which would adversely impact the change detection output leading to a misleading interpretation of the phenomenon (deformation) that one seeks to measure and observe.
Consider the visual output of this step the same as Figure 5 as there is no discernible change to it after running the Apply Orbit Files operator. .
d) Back Geocoding - Now that I have updated the Orbit files, the Back Geocoding operator will help me to assign the energy reflectance values captured by the satellite's receiver to each pixel representing the Earth's surface. This conversion of data, originally acquired in radar coordinates to regular coordinates necessitates the use of a Digital Elevation Model. DEM is a dataset of the elevation values of the bare Earth surface - these exclude the elevation data of natural features such as vegetation or built features such as buildings which lie above the earth's bare surface.
The Back Geocoding operator not only assigns the reflectance values to pixels across both the Imagery datasets but also it bundles both the outputs in a single Stack (as seen in the top-left section of Figure 6 below). This bundling step is technically known as Coregistration and has a separate operator by the same name within SNAP if one wants to perform just this operation alone.
Coregistering SAR Imagery datasets is a mandatory step as only upon doing so will both the datasets become comparable for SNAP software i.e. suitable for joint analysis for change detection purposes, something which I also need to do for this Deformation Mapping workflow.
e) Enhanced Spectral Diversity - While the previous operation - Back Geocoding - involved interpolating the data acquired in radar geometry onto the surface of the earth, this operator utilizes an optimization algorithm to improve the geometric precision of the assignment of energy reflectance values onto pixels in the SAR Imagery.
This is the last step in the 'Coregistration' phase of the workflow. The Coregistered ESD Stack now contains totally comparable sets of Imagery pre and post eruption - ripe for applying Interferometry.
f) Interferogram - This operator is the first step in my quest to detect Deformation and measure the Ground Displacement rate between the pre and post eruption SAR Imagery datasets, with centimeter-level accuracy.
The three Interferometric outputs generated by this operator are the Intensity, Phase and Coherence bands-
Intensity - By default, as applicable to the Intensity band of an individual Imagery dataset (Figure 4), the Intensity values are a square of the Amplitude values where Amplitude is the reflected microwave energy from the Earth's surface, originally transmitted by and subsequently captured by the satellite itself. However, in the Interferogram output for the Coregistered Stack, the Intensity band is formed in a different way - by the multiplication of the Amplitude value for the same pixels across both the datasets i.e. Intensity of Pixel 1 in Interferogram output of Coregistered Stack = Amplitude of Pixel 1 in Imagery 1 x Amplitude of Pixel 1 in Imagery 2.
Phase - Likewise, by default and as is applicable the Phase band of an individual Imagery dataset (Figure 4), Phase value of a pixel is a modified representation of the distance between the satellite's antenna and the ground target. As I am seeking to measure Displacement (rate of uplift or subsidence), the Phase value or rather, the change in Phase value across both the datasets is what matters. The Phase band in the Interferogram output is exactly this - it captures the Phase difference i.e. derived by subtracting the Phase value of a pixel in one dataset from the Phase value of the same pixel in the other dataset - Phase of Pixel 1 in Interferogram output of Coregistered Stack = Phase of Pixel 1 in Imagery 1 - Phase of Pixel 1 in Imagery 2. The Phase output is the 'Interferogram'.
Coherence - Coherence is a measure of Interferogram quality and the Coherence of a pixel is a normalized value ranging from 0 to 1 which is derived from the default values of Amplitude and Phase across both the individual Imagery datasets. If the Amplitude value for a pixel across both the Imagery datasets are similar and the Phase value for the same pixel across both the Imagery datasets are similar too, then the Coherence computation takes a value closer to 1 i.e. Very High Coherence. If the Amplitude and/or Phase values across both the datasets are dissimilar, then the Coherence value diminishes towards 0 i.e. Low Coherence. What a Deformation researcher is keen to observe in the Interferogram output of the Coregistered Stack is the latter - clusters of pixels with Low Coherence, which is indicative of considerable difference in Amplitude and/or Phase values across both the datasets, a telltale sign of Displacement.
g) TOPSAR Deburst - TOPSAR mode of Imaging, as used by Sentinel-1, involves the transmission of microwave energy by a moving satellite in short bursts which illuminates sections of the ground target (Figure 4). The data discontinuity between two bursts result in the formation of horizontal black lines between the Imagery sections / sub-swaths. The Deburst operator fills in the missing data by estimating the values (technique described here), thereby making the entire Imagery seamless.
h) Topographic Phase Removal - The Phase in the Interferogram output of the Coregistered stack (i.e. the difference in Phase between both the individual Imagery datasets), occurs due to two reasons-
1) change in Topographic features (eg. flooding, harvesting, wind reorienting vegetation etc.), and 2) movement (uplift / subsidence) of ground surface i.e. Displacement
In order to visualize Deformation, it becomes important to remove the former - the Topographic Phase, so that what is left is just the Phase (difference) due to Displacement.
i) Multilooking - This operator utilizes an algorithm to reduce the inherent Speckle Noise in the Interferogram output by standardizing all the pixels, including the ones with inconsistent sizes caused by the angle of the incoming energy signal, into similar-sized 'square' shapes. Thus, Multilooking helps to improve the radiometric resolution of the output at the expense of degrading its spatial resolution. The clarity of the Interferogram output improves from the application of this operator, as evident in Figure 11 below, thereby allowing for easier loss-of-Phase detection and interpretation.
j) Goldstein Phase Filtering - This and the next step is where it gets too technical for me to both, understand as well as explain what is happening in a simple paragraph or two. To summarize from what I have gathered - like Multilooking, Phase Filtering is also used to reduce the noise in the Interferogram output (compare Figure 12 with Figure 11). Using a Phase Filter right now also helps in generating a more accurate output in the next step involving Phase Unwrapping.
k) Snaphu Export - Snaphu (Phase Unwrapping) - Snaphu Import - Phase to Displacement - Technically, Snaphu (Statistical-cost, Network-flow Algorithm for Phase Unwrapping) is a critical step in the Interferometric processing chain - it involve the execution of a complex algorithm.
While the cleaned-up Interferometric Phase gave us an indication of where the Deformation was occuring (by observing the location of fringes i.e.), Phase Unwrapping using Snaphu and subsequently, the 'Phase to Displacement' operator helps in the calculation of precise Displacement rates between the two sets of Imagery pre and post the eruption event.
l) Range-Doppler Terrain Correction - Besides removing the Displacement values over non-land features through the use of Digital Elevation Model (DEM), the primary objective of running the Range-Doppler Terrain Correction operator is to make the Interferogram output geometrically accurate (orthorectification) to the Earth's globular shape. The tilt of satellite, the tilt of its sensor, curvature of the earth among other aspects distort the spatial relationship between the pixels and Terrain Correction helps the output to resemble the real world i.e. how it would appear on a globe.
m) Write (& Export) - This operator enables me to download the output on the system for further processing / visualization purposes.
Visualizing the Interferometric Phase & Displacement Output using Google Earth
Phase Visualization on Google Earth
Video 2: Fogo (Cape Verde) Volcano's Interferometric Phase layer has been overlaid on Google Earth Basemap. Strong Deformations (Fringes) are present around the cone of the volcano - both Subsidence as well as Uplift areas can be observed (refer Color Scale)
Video 3: Mt. Nyiragongo (Congo) Volcano's Interferometric Phase layer has been overlaid on Google Earth Basemap. Strong Deformations (Fringes) can be observed - however, these are not around the cone of the volcano, but farther away.
What could these Deformations be attributed to, then?
These are due to persistent seismic activity (Earthquakes) following the eruption.
Displacement Visualization on Google Earth
Video 4: From Fogo (Cape Verde) Volcano's Displacement Map, one can clearly observe a) the location of Deformations, b) the type of Deformation / Displacement - Uplift or Subsidence (refer Color Scale) and c) the quantitative rates of Displacement.
Thank you for staying on till the end of this comprehensive article. Do watch the Video walkthrough in case you'd like to observe the way to use SNAP software to execute the processing chain to map Deformation at Volcanic sites.
Finally, let me leave you with a spooky possibility so that you do not crave for more volcanoes😁.
ABOUT US
Intelloc Mapping Services | Mapmyops.com is based in Kolkata, India and engages in providing Mapping solutions that can be integrated with Operations Planning, Design and Audit workflows. These include but are not limited to - Drone Services, Subsurface Mapping Services, Location Analytics & App Development, Supply Chain Services, Remote Sensing Services and Wastewater Treatment. The services can be rendered pan-India, some even globally, and will aid an organization to meet its stated objectives especially pertaining to Operational Excellence, Cost Reduction, Sustainability and Growth.
Broadly, our area of expertise can be split into two categories - Geographic Mapping and Operations Mapping. The Infographic below highlights our capabilities.
Our 'Mapping for Operations'-themed workflow demonstrations can be accessed from the firm's Website / YouTube Channel and an overview can be obtained from this flyer. Happy to address queries and respond to documented requirements. Custom Demonstration, Training & Trials are facilitated only on a paid-basis. Looking forward to being of service.
Regards,