INSAR-BASED LOS SURFACE DEFORMATION COMPARISON OF METRO MANILA BEFORE AND AFTER THE JANUARY 2020 TAAL VOLCANO ERUPTION

The city of Metro Manila has been constantly battered by several hazards on an annual basis. On January 2020, the Taal Volcano erupted with multiple recorded earthquakes. Previous literatures have found that Metro Manila is experiencing a steady subsidence. Determination of land uplift or subsidence is crucial in planning and mitigating the effects of flooding in the area. The study aims to determine whether an uplift occurred in Metro Manila after the eruption or is the study area still experiencing subsidence This study uses a pair L1 SLC Sentinel 1 Images. Radar Interferometry is used to generate Interferograms and Satellite Line of Sight (LOS) deformation was determined between the 2 dates of image acquisition. It was found that the Metro Manila area generally experienced an uplift except for some areas in Caloocan which shows subsidence. The uplift magnitude gradually decreases going from the South to North with a max value of 9.6cm.


BACKGROUND
Every year, the Philippines experiences an average of 20 typhoons (Padagdag, 2018). Compounding with climate change, the magnitude of rainfall brought about by these typhoons may result in extreme flooding. The flooding events cause enormous casualty, economic and environmental losses. Metro Manila is one of the key areas where these flooding events takes most of the toll. It is suggested that 3 decades from now, Metro Manila will be submerged due to the effects of these extreme weather events (Ng, 2020).
One of the main risk factors when creating flood mitigation program is the determination of accurate flood models. This can be done by generating an accurate elevation model of the study area. The challenge occurs due to the dynamic nature of the topography. Most especially due to the presence of Surface Deformationparticularly, Subsidence and Uplift (Eco, 2020). Subsidence and Uplift are defined as the loss and increase of surface elevation, respectively. These change in elevation can come from a variety of sources. In an article by Lagmay (2011), Metro Manila is generally experiencing a subsidence due to groundwater extraction.
Without accounting for this deformation, the flood mitigation programs will gradually become less effective and will need to be updated every so often. This ensures that the change in the topography is accounted for when performing the flood modelling to generate an accurate risk assessment that would translate to better flood mitigation techniques and programs.
Last January 12, 2020, the Philippine Institute of Vulcanology and Seismology (PHIVOLCS) began showing signs of unrest. A total of 673 volcanic earthquakes were recoded within the timeframe (UN Migration, 2020). Phreatic explosions with Ash Plums were also observed and alert level 4 was raised due to the continuous eruptions. The province of Batangas was placed under a State of Calamity 2 days after. A total of 61,123 households were affected and 235,655 individuals were displaced and evacuated.
The advent of remote sensing brought the technology of measuring the earth's surface without having to do actual ground fieldwork. Synthetic Aperture Radar (SAR) is an active remote sensor. It has been used to map out deformation changes on the earth surface. Radar satellites are good with measuring distances. These distance measuring capabilities are what makes Radar Remote Sensing a useful tool in detecting Surface Deformation (Woodhouse, 2006).
In this study, pairs of Sentinel 1-A L1 SLC Images will be used. Line of Sight (LOS) Deformation changes will be measured between image pairs and will be measured relatively. The concept of Interferometry will be used to determine if an uplift or subsidence occurs in the area.

OBJECTIVES
The objective of this study is quantifying the surface deformation experienced by Metro Manila due to the Taal Volcano Eruption. Previous literature has analysed multiple time series SAR data and found that Metro Manila generally experience subsidence (Eco, 2020). The study aims to see if the The specific objectives are the following: 1. Determine presence of subsidence or uplift in the study area after the eruption 2. Quantify amount of LOS Deformation in the study area 3. Determine which areas are affected the most in terms of the LOS Deformation

LITERATURE REVIEW
Some of the most common and reliable methods of measuring surface deformation or elevation changes are the use of Field Levelling, GNSS and DInSAR (Sneed, n.d.). These methods have their own advantages and tradeoffs depending on the accuracy and extents of the measurement. Factors such as costs, timeframe, manpower and availability of the equipments play heavily on which methods should be used.
SAR is an active remote sensing technique. It uses an antenna that emits an electromagnetic wave in the microwave region that interacts with the objects on the ground then gets backscattered to the antenna (Meyer, n.d.). The returning signal determines the distance of the object from the antenna and introduces some additional information about the feature on the ground as a function of the amount of energy backscattered to the antenna (Liew, 2001). The strong capability of SAR sensors to determine distance measurements is what makes the technology viable for determining surface deformations. This technique is suitable for large scale deformation due to the vast amount of data it can gather at a single pass (Sneed, n.d.).
SAR Interferometry (InSAR) is a technique that looks at the phase information between 2 or more SAR Images (Goldstein, 1989). The phase difference between these images is interpreted as the sum of all the phase contributions in the two images. Alternatively, it can also be viewed as the change in distance between the two time frames where the SAR images were acquired (Massonet, 1998).
InSAR works by deriving a distance measurement from Image 1 then detects another distance measurement from Image 2. In Figure 1, the distances Z from the Satellite 1 and 2 were derived from different time periods. This is more specifically known as a Repeat Pass Interferometry wherein the satellite passes in the same area after a period of time (Matsuoka and Yamazaki, 2000). In the case of Sentinel 1, it passes by the same area after 12 days.
The two distances are compared by generating the Phase Difference between the 2 signals. The Phase Difference is represented by an Interferogram which can be obtained by multiplying the phase of Image 1 with the complex conjugate of Image 2. The resulting Interferogram consists of all the phase contributions of different sources in the SAR Images including the phase difference due to the deformation. It is important to note that the values in the Interferogram are wrapped from 0 to 2π or one cycle of the wavelength used by the specific sensor (Crosetto, 2002). The quality of the phase difference measurement is dependent on whether drastic changes may have occurred in the area during the timeframe of the data acquisition. To quantify this, a measurement called the Coherence is derived. This variable describes the statistical similarity of a pixel taken from Image 1 and Image 2. The value of the Coherence will determine whether the deformation contribution within the Interferogram can be retrieved properly. Otherwise, the signal is plagued with noise and deformation information cannot be reliably measured (Zebker, 1994). Once the coherence is resolved, unwrapping of the Interferometric Phase in the Interferogram is done to derive the actual linear measurement of the deformation (Gudmundsson, 2002). Though, large in scale in terms of the data gathering, the accuracy derived from a pair of SAR images in terms of deformation measurement are not at par with the likes of GNSS and Field Levelling (Serrano-Juan, 2017) where levelling still remains to be the standard for vertical measurements. However, several studies have already shown the effectiveness of using SAR in geohazard studies (Tomas & Li, 2017) such as landslides, tectonic movements (Pepe & Calo, 2017) and deformations due to volcanic activities (Zhou, 2009).

Data Pre-processing
Two Sentinel 1A SLC Images taken 12 days apart are utilized in this study. Images 1 and 2 were taken last January 9, 2020 and January 21, 2020, before and after the Jan.12 Taal eruption. Level 1 SLC products are chosen since these datasets retain both the Amplitude and Phase component of the returning signal.
Application of Precise Orbit (POE) files to the images are done. This is so that the stacking of the images can be done accurately and a pixel-to-pixel comparison between both images can be done. TOPS Debursting is then applied to select the specific burst/s that contain the study area. TOPS Split is used to further crop the stack to the study site only.

Interferogram Generation
An Interferogram represents the phase difference between the 2 SAR images. This is computed on a pixel-to-pixel basis by multiplying the phase of Image 1 with the complex conjugate of the phase of Image 2. This is calculated in a per pixel basis given by equation 1  Interferogram values are wrapped within 0 to 2π since only cyclic difference between the phases are measured in the Interferogram. Each color cycle represents a deformation value corresponding to 1 wavelength. The spatial variation of these phase differences is presented as fringes.
(1) The resulting Interferogram is a sum of the different phase contributions from difference sources. The first three contributors are Flat Earth Phase, Topographic Phase and Deformation Phase. These phase contributions are deterministic in nature. Removal of the Topographic and Flat Earth Phases results in a Differential Interferogram. The next contributor is the atmospheric phase. The timeframe considered for this study will be sufficient to assume that there has been no significant atmospheric change between the two images. The last contribution is the noise phase or the Speckle.  This is directly related to Coherence. The higher the coherence, the lower the contributions of the noise phase in the Interferogram.

Phase Unwrapping
Interferogram values are wrapped within 0 to 2π since only cyclic difference between the phases are measured in the Interferogram. Each color cycle represents a deformation value corresponding to 1 wavelength. To get the actual linear measurement, the process of Unwrapping is done. In this study, the Minimum Cost Flow (MCF) unwrapping algorithm is utilized. This is done by first selecting a seed point where the phase difference relative to the seed point is expanded and absolute phase difference is derived for all succeeding points. This process allows for the determination of the actual phase difference accounting for the cycle ambiguity to be counted. Phase Unwrapping of the image stack is done through the Statistical-Cost, Network-Flow Algorithm for Phase Unwrapping (SNAPHU). After Unwrapping is done, Terrain Correction is performed (Serco Italia SPA, 2018). the resulting image will be overlain with the boundary of the study area to visualize the results.

RESULTS
To reduce the processing time, the images in figure 2 were split using the TOPS Split function. This allows the removal of unwanted bursts in the image that will no longer be part of the processing. Figure 6 are the results after splitting. .

Figure 6. TOPS Split SAR Images
Six (6) bursts were removed from the original image. These two resulting images are coregistered using S-1 Back Geocoding with an SRTM 30m DEM. After coregistration, ESD is applied and the interferogram is generated. See figure 7 for the interferogram.

Figure 7. Interferogram
The interferogram is then filtered using Multilooking to reduce speckle and create square pixels. An adaptive filter, Goldstein Phase Filter, is applied. This filter that boosts local dominant fringe patterns. The results of these filters are seen in figure 8. On the opposite end, highly urbanized areas are suitable for interferograms because they have a high coherence value. This The International Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences, Volume XLVI-4/ W6-2021Philippine Geomatics Symposium 2021, 17-19 November 2021 is due to the fact that signals backscatter strongly in urban areas due to the strong presence of double bouncing and that urban areas are highly unlikely to experience any sudden change that may affect the backscattered signal during the succeeding image acquisition. The greater the number of color cycle changes, the greater the surface deformation that may have occurred.
Terrain correction is applied to the filtered interferogram so that the image will now have a map geometry and proper coordinate reference system. Figure 9 shows the terrain corrected interferogram overlain in Google Earth.

Figure 9. Terrain Corrected Interferogram
Unwrapping of the interferogram is done to obtain the absolute phase difference using the SNAP plugin SNAPHU. This is a plugin that allows for phase unwrapping using the MCF algorithm. Figure 10 shows the unwrapped phase of the study area overlain in Google Earth.

Figure 10. Unwrapped Interferogram
Unwrapping of the interferogram is done on a per tile basis. This variable is user specified. For this study, 10 tiles per row and column were the default. The tile is unwrapped one at a time until the entire image is processed. After unwrapping is done, the resulting image is imported back to SNAP and can now be converted to displacement. Note that the displacements measured in the study is in terms of the satellite Line of Sight (LOS). Interferometry measures the change in distance relative to the line of sight of the sensor. Figure 11 is the converted unwrapped phase to displacement measurements.   Figure 13 shows that a month prior to the eruption, the city of Metro Manila is experiencing subsidence centered around San Juan City with portions of Manila, Quezon City, Makati and Mandaluyong and is spreading radially outward. This image shows that prior to the eruption, Metro Manila is subsiding relative to the SAR LOS with values of up to -2.5cm. This somewhat confirms the previous studies that, generally, Metro Manila is subsiding in terms of surface deformation.
A month after the eruption, as seen in figure 14, Metro Manila is, once again, experiencing subsidence with a general behaviour much like in Figure 12 but with opposite directions. Areas that experience the highest uplift in the original data now experienced subsidence where the magnitude gradually decreases towards the north. Areas with the greatest subsidence include Muntinlupa, Las Pinas and Paranaque with a magnitude of -5.00cm to -7.5cm. Areas within the central Metro Manila are also experiencing subsidence with magnitude of around -2.5cm with no areas experiencing any uplift.
Both figures 13 and 14 have been filtered based on coherence values. The higher the coherence, the more reliable the deformation measurements are. Unlike figure 12 which covers the entire area of Metro Manila, figures 13 and 14 have some areas particularly in the upper right portion of Quezon City and Caloocan which covers the eco park mentioned previously. The coherence values within those areas are not enough to estimate a proper deformation measurement. This area also coincides with the areas in figure 12 where subsidence is observed. Thus those measurements should not be interpreted directly due to the low coherence value in that area.

CONCLUSIONS AND RECOMMENDATIONS
The study was able to show the presence of surface deformation in the area using the technique of InSAR. The study was also able to show how InSAR can be used to collect large amounts of data without having the need for an actual fieldwork. In terms of resources, moderate computational capacity is needed to perform the analysis done in the studywhich makes the technology all the more viable for monitoring purposes as new Sentinel images are generated every 12 days. Factoring in both the ascending and descending orbits of 1A and 1B, the temporal resolution can be reduced to 6 days only.
The maximum uplift and subsidence detected was 9.6cm and -4.7cm respectively. Generally, the study area experienced an uplift with a little portion experiencing some subsidence. The amount of surface deformation was quantified in terms of the subsidence and uplift that were present in during the time frame of the data acquisition. The areas where the uplift and subsidence were also determined by mapping out the extents of the surface deformation.
The study only made use of 1 pair of SAR images. More time series analysis is preferred to get a better understanding of the behavior of the surface deformation in the area.
For this study, it is recommended to use multiple SAR images for create a larger stack for better analysis. Furthermore, the technique of Persistent Scatter Interferometry (PS-InSAR) will greatly benefit this study. Since the study area is mostly comprised of urbanized areas, possible sudden changes that may greatly affect the Coherence of the signal will be minimized thus can result in a better analysis of PS points when using PS-InSAR.