Automated Methods for Detecting Volcanic Deformation Using Sentinel‐1 InSAR Time Series Illustrated by the 2017–2018 Unrest at Agung, Indonesia

Radar satellites, such as Sentinel‐1, are now able to produce time series of ground deformation at any volcano around the world, but atmospheric effects still limit the real‐time detection of unrest at tropical volcanoes. Here, we test two approaches to correct atmospheric errors—phase elevation correlations and global weather models—and assess the ability of Interferometric Synthetic Aperture Radar (InSAR) time series to detect deformation anomalies using either a fixed threshold or a cumulative sum control chart. We use the 2017–2018 crisis at Agung volcano as a case example because strong atmospheric signals were originally misidentified as true deformation, and obscured the subtle deformation pattern associated with magmatic activity. We assess the Receiver Operating Characteristics (ROC) of each method and found the average area under the ROC curve to be about 0.5 for the uncorrected data (corresponding to no discrimination capability), around 0.8 after combined atmospheric corrections (weather model and phase elevation approaches), and more than 0.95 using a cumulative sum control chart (where 1 corresponds to ideal separation between classes). Our results retrospectively show that uplift could have been detected to a 95% level of confidence for both ascending and descending time series by October 2017, 15 days after the start of the seismic swarm and 1 month prior to the eruption. Thus, our approach successfully flags anomalous behavior without relying on visual inspection or selection of an arbitrary threshold, and hence shows potential as a monitoring tool for volcano observatories globally.


Introduction
Radar interferometry (InSAR, Interferometric Synthetic Aperture Radar) has been used successfully to quantify surface ground deformation on volcanoes located in a range of tectonic settings including subduction arcs (Lu & Dzurisin, 2014;Pritchard & Simons, 2002), continental rifts (Biggs et al., 2009;Wright et al., 2006), and oceanic basaltic islands (Froger et al., 2004;Lundgren et al., 2013;Sigmundsson et al., 2015). However, real-time volcano monitoring largely relies on ground-based measurements of deformation (e.g., Global Positioning System, GPS, or tiltmeters) (Fiske & Kinoshita, 1969;Linde et al., 1993;Larson et al., 2010;Voight et al., 1998), and satellite InSAR has typically only been used retrospectively. However, the latest generation of InSAR satellites have much shorter repeat times than previously available, meaning that images are now collected sufficiently frequently to be useful in volcanic crises (Ebmeier et al., 2018). Recent examples include the regional study of volcanic unrest in Latin America (Pritchard et al., 2018) and the study of specific events such as the 2014 unrest at Chiles-Cerro Negro  or the 2010 eruption at Merapi volcano (Pallister et al., 2013). Automated monitoring systems are being developed to routinely analyze large volumes of InSAR data, making use of artificial intelligence tools such as machine learning (Anantrasirichai et al., 2018(Anantrasirichai et al., , 2019 and blind source separation (Ebmeier, 2016;Gaddes et al., 2018Gaddes et al., , 2019. However, the accuracy of InSAR measurements is highly variable from one volcanic environment to another (Ebmeier et al., 2013;Parker et al., 2015) and this places a limitation on the use of InSAR as a monitoring tool. Deformation measurements are affected by several factors: (1) temporal decorrelation induced by changes of surface scattering (e.g., water or vegetation), (2) SAR geometry, which produces layover and shadowing in high relief areas, and (3) atmospheric artifacts, which include interactions with the ionosphere and temperature/pressure variations in the troposphere, but are dominated by variations in tropospheric water vapor content. Andesitic stratovolcanoes located on tropical islands are especially challenging environments for InSAR measurements as they combine all these limiting factors: vegetation cover, steep topography, and In this paper, we present an innovative strategy to correct atmospheric delays, specifically targeted at steep, tropical volcanoes. The method combines high-resolution weather models with a phase elevation approach to reduce the atmospheric signals. We produce corrected InSAR time series that can be used to automatically detect sequential anomalies related to ground deformation signals. We test the various approaches available using the 2017 crisis at Agung volcano (Bali, Indonesia) as a case example and find that detection of anomalous deformation would have been possible 1 month prior to the eruption. Our methodology is suitable for implementation in automated InSAR Sentinel-1 processing systems, which will benefit volcano observatories who are routinely monitoring volcanic unrest.

Empirical Approach
Tropospheric phase delays occur in interferograms due to changes in water vapor, pressure, and temperature in the troposphere between acquisitions. The delays are typically decomposed into a stratified component, which is related to vertical stratification of atmospheric properties and hence correlates with the topography and a turbulent component, which is spatially correlated over short wavelengths and temporally uncorrelated (Hanssen, 2001). Correction approaches can be divided into three categories: those based on empirical relationships (e.g., phase elevation), those based on external data sets (e.g., GPS, MODIS), and those based on weather models (e.g., ECMWF).
The simplest approach to correcting stratified delays is to estimate the correlation between the phase of the interferogram and the elevation given by a digital elevation model (DEM) assuming linear correlation (Beauducel et al., 2000;Remy et al., 2003) or power law model (Bekaert et al., 2015). The phase elevation correlation can be calculated either over the entire interferogram (Cavalié et al., 2008) or on local windows to take into account the spatial variability of the signals (Béjar-Pizarro et al., 2013).
Although the method is easy to implement and provides good results at the first order, its application is limited in volcanic regions as deformation signals often correlate with the topography (Delacourt et al., 1998;Ebmeier et al., 2013). Furthermore, this approach does not correct the turbulent component of the atmospheric signal, which can be significant in tropical regions.

External Data Sets
To address the limitations of the empirical method, external data sets have been used to calculate troposheric phase delay maps to correct atmospheric signals. Data sets include (1) local meteorological data (Delacourt et al., 1998), (2) spectrometer data acquired by the MODerate resolution Imaging Spectroradiometer (MODIS) (Li et al., 2005(Li et al., , 2009 or the MEdium Resolution Imaging Spectrometer (MERIS) (Li et al., , 2012, and (3) zenith tropospheric delays (ZTD) estimates from Global Navigation Satellite System (GNSS) sites Onn & Zebker, 2006;Williams et al., 1998).
The MODIS spectrometer measures water vapor delay maps at high spatial resolution (0.3-1.2 km) with good accuracy (root-mean-square of ∼1 mm compared to GPS) (Li et al., 2009), but this approach only works in cloud-free conditions and the time of the acquisition rarely coincides with the time acquisition of Sentinel-1 and ALOS-2 data (>5-hr difference).
GPS tropospheric delays have a better temporal resolution, but the spatial resolution is highly dependent on the density of the GNSS network. Although the benefit of such methods has been demonstrated in places with dense network such as California or United Kingdom (Yu et al., 2017, its application to the study of volcanoes is limited as many tropical volcanoes are instrumented poorly or not at all.

Journal of Geophysical Research: Solid Earth
10.1029/2019JB017908 or U.K. Met Office's Unified Model have been successfully applied to derive high-resolution (<10 km) atmospheric delay maps over volcanic area such as Hawaii or Etna (Foster et al., 2006;Wadge et al., 2010). Nico et al. (2011) have also simulated 1-km resolution atmospheric phase delays using the Weather Research and Forecasting model. However, such high-resolution local weather models are only available for limited geographical areas.
For this reason, many InSAR studies have preferred global weather models based on the reanalysis of meteorological data, including the Re-Analysis ERA-Interim from the European Centre for Medium-Range Weather Forecasts (ECMWF) (Doin et al., 2009;Jolivet et al., 2011) or the North American Regional Reanalysis from the National Centers for Environmental Prediction/National Center for Atmospheric Research (NCEP/NCAR) (Pinel et al., 2011;Parker et al., 2015).
ERA-Interim corrections have been successfully applied to correct long wavelength atmospheric signals on tectonic areas (Jolivet et al., 2011), but the coarse spatial resolution (∼75 km) is less suitable for correcting atmospheric signals over the smaller scale of individual volcanic edifice (∼1-30 km) (Grosse et al., 2009) than higher-resolution models such as the North American Regional Reanalysis (∼32 km) (Parker et al., 2015). Moreover, the release of ERA-Interim products is delayed by several months and therefore they cannot be integrated into a real-time InSAR processing strategy for volcano monitoring. ERA5 has now replaced the ERA-Interim reanalysis, which stopped being produced on 31 August 2019. The ERA5 products will have higher temporal (1 hr) and spatial (∼30 km) resolution than ERA-Interim and preliminary estimates would be available within 5 days. Therefore, ERA5 products could be tested in the future for performing atmospheric corrections over volcanic regions.
In the meantime, a high-resolution ECMWF global forecast (HRES-ECMWF) was released in 2016 with a spatial resolution of 9-12 km and a temporal resolution of 6 hr. For example, such products have been successfully used in the northern Tibet for correcting Sentinel-1 interferograms from strong tropospheric signals to extract coseismic and interseismic displacements (Shen et al., 2019). In addition, a Generic Atmospheric Correction Online Service (GACOS) applies an Iterative Tropospheric Decomposition (ITD) to the HRES-ECMWF products to separate and interpolate the stratified and the turbulent components of the tropospheric delay , Yu, Li, Penna & Crippa 2018. Thus, GACOS provides high-resolution phase delay maps only 1-2 days after the Sentinel-1 acquisition, which are suitable for producing corrected InSAR time series for monitoring purposes. The GACOS corrections perform well in test case studies from the United Kingdom, California, New Zealand, Italy, Tibet, Nepal, Iceland, and Algeria , but their performance in the more challenging case of a tropical volcanic environment has yet to be assessed. Here, we test the suitability of GACOS atmospheric corrections for near-real-time deformation monitoring with the aim of delivering corrected time series of ground deformation less than 2 days after the acquisition. The method is tested at tropical volcanoes using the example of the 2017-2018 volcanic crisis at Agung volcano, Indonesia .

The 2017 Crisis at Agung Volcano, Bali
Agung is a stratovolcano (3,031 m above sea level) located on the island of Bali in Indonesia (Figure 1a). After more than 50 years of inactivity since a VEI5 eruption in 1963 (Self & Rampino, 2012), the volcano showed signs of reawakening between September and November 2017 . The unrest included an earthquake swarm recorded by seismic stations operated by the Indonesian Center for Volcanology and Geological Hazard Mitigation, which was felt by the local population as well as visual observation of degassing. The unrest lasted almost 3 months and culminated in an eruption on 21 November 2017. During the 1963 eruption, more than 1,000 people were killed (Witham, 2005), whereas during the 2017 crisis, ∼140,000 people were evacuated, preventing any fatalities Syahbana et al. (2019).
The ability to provide scientific advice during volcanic crises is highly dependent on the quality of the monitoring data and the ability to identify the underlying physical processes or forecast behavior. Agung has only had two eruptions in the last century and is not as well monitored as more active Indonesian volcanoes, such as Merapi (Pallister et al., 2013). At the time of the 2017 seismic crisis, the ground network at Agung consisted of four seismic stations and five GNSS stations, which were insufficient to provide precise information about the geometry, location, or evolution of the magmatic sources.
Public communication is critical during volcanic crisis and social media (e.g., Twitter, Facebook, WhatsApp) were intensively used during Agung crisis to disseminate information to the population (Syahbana et al.,  Figure 1b). However, a GACOS estimate for the atmospheric delay on this interferogram visually demonstrated that the fringes were caused by an atmospheric artifact rather than deformation ( Figure 1c). The misinterpretation of atmospheric artifacts as volcanic ground deformation is not unusual in InSAR analysis and has been previously reported at Etna and Llaima volcanoes (Beauducel et al., 2000;Remy et al., 2015). The diffusion in real time of inaccurate information had a large negative impact on the management of the Agung volcanic crisis (e.g., self-evacuation or people selling their livestock at a loss)  and highlighted the potential dangers for volcano observatories of open-access data policies and social media.
During the crisis, we use GACOS to correct individual 12-day interferograms as they became available and the residual on Agung edifice following correction was usually small (less than one fringe). Thus any deformation was initially considered to be below the detection threshold. By November, a persistent low-magnitude residual signal had been noted in descending track data on the NE flank, but as no corresponding signal was seen in the ascending track data, this was dismissed as recurrent atmospheric turbulence associated with the coastal setting. Retrospectively, atmospherically corrected InSAR images were used to detect a broad, preeruptive deformation signal located on the NE flank of Agung volcano, which has been interpreted as a dike intrusion connecting the edifices of Agung and Batur . In this paper, we assess whether near-real-time atmospheric corrections to time series data could have made the Sentinel-1 data more useful for managing the crisis.

Sentinel-1 Data
Sentinel-1 is a radar mission operating at C-band ( = 5.6 cm) launched in 2014 by the European Space Agency. Using the default TOPS mode, radar images are acquired with a swath width of 250 km, enabling large-scale regional mapping (Geudtner et al., 2014). At the time of writing, the constellation is composed of two satellites, Sentinel 1A and Sentinel 1B, which together have a repeat cycle of 6 days over Europe and 12 days for the rest of the world. In terms of volcano monitoring, this is a significant improvement over the irregular acquisition strategy and 35-day repeat interval of previous C-band missions such as ERS-1/2 and ENVISAT. Due to its system capabilities and policy of freely available data, the Sentinel-1 mission opens new possibilities for routine monitoring in a range of fields, including the monitoring of ground deformation associated with natural hazards (e.g., earthquakes, volcanoes, or landslides) or anthropogenic activities (e.g., water pumping or oil/gas extraction). However, as with previous C-band sensors, short time span interferograms generated with Sentinel-1 data will be dominated by tropospheric signals that must be accounted for when designing monitoring algorithms.
We processed Sentinel-1 SAR images of Bali ( Figure 1a) using GAMMA software and produced a total of 31 ascending interferograms and 31 descending interferograms for a 1-year period (from May 2017 to May 2018). This includes ∼6 months prior to the eruption on 21 November 2017 and ∼6 months after. As expected in tropical regions, strong atmospheric signals are present in the majority of the 12-day wrapped interferograms, which mask the detection of ground deformation (Figures S1a and S2a in the supporting information).

Atmospheric Corrections
We test two near-real-time approaches for correcting atmospheric signals: (1) empirical phase elevation corrections (Beauducel et al., 2000;Bekaert et al., 2015;Remy et al., 2003) and (2) the latest generation weather models using the GACOS , Yu, Li, Penna & Crippa 2018. In this context, we will be able to compare the two approaches and to quantify their relative performance.
For the simple empirical approach, we calculate the linear correlation between the phase of each interferogram and the elevation of the 90-m DEM provided by the Shuttle Radar Topographic Mission. The phase elevation correlation is estimated on the entire interferogram rather than on the edifice alone to reduce the chances of removing small ground deformation signals that may correlated locally with the edifice topography.
GACOS provides Zenith Tropospheric Delay (ZTD) maps based on the HRES-ECMWF weather model (0.125 • spatial resolution and 6-hr temporal resolution) and interpolated to a 90-m grid using the Shuttle Radar Topographic Mission DEM and the ITD algorithm . ITD can be applied on the pointwise GPS ZTD estimates, the regularly spaced HRES-ECMWF grids, or a combination of both. Compared with the classic DEM mapping of the tropospheric delay, the iteration procedure improves the model of elevation-dependent signals, especially over mountainous areas (Yu et al., 2017) and during rainy seasons ). The GACOS system can be accessed by a command-based API, which is fully automatic without human interactions. Therefore, it is possible to routinely acquire the ZTD maps without having to download, decode, and interpolate the HRES-ECMWF.
We calculate phase difference between two GACOS ZTD maps obtained at the time of Sentinel-1 acquisitions and project the result into the corresponding line-of-sight (LOS) using the incidence angle of the SAR image. Because GACOS atmospheric delays and interferograms are relative measurements, we chose a common reference point far from the volcanic center. After the GACOS corrections, we notice the presence of a residual phase ramp in a few interferograms, which we remove using a best fit plane. Figure 2a shows the 12-day descending track interferogram from 28 August 2017 to 9 September 2017 (20170828-20170909). The initial standard deviation is about 2 cm, and the atmospheric signal consists of a long wavelength positive LOS change located along the N-NE coast of Bali. In this case, the GACOS weather model explains the observed pattern well (Figure 2b). Although part of the signal correlates at first order with the topography (Figure 2d), the linear phase elevation approach is not able to correct the atmospheric delay ( Figure 2c). The standard deviation after the GACOS correction is 1.2 cm (Figure 2e), which is better than the 1.8 cm obtained with the phase elevation approach (Figure 2f). After the GACOS corrections, there is a persistent negative signal at Agung volcano at elevations above 1,700 m (Figures 2e and 2d). Because Agung (3,031 m) is the highest peak on Bali (Figure 1a), we expect local atmospheric perturbations (wavelength less than 10 km) with different characteristics to those over the rest of the island, which could explain why both approaches fail to correct such signals. Figure 3a shows the 12-day ascending track interferogram from 16 to 28 May 2018 (20180516-20180528). The initial standard deviation is 1.1 cm, and the atmospheric signal consists primarily of small positive LOS changes around the volcanic peaks, Batur and Agung. Although the GACOS weather model predicts the regional pattern of the atmospheric signal well (Figure 3b), the correction introduces some large residuals as it overestimates the phase delay in the NW corner of the interferogram (Figure 3e). The standard deviation increases to 1.7 cm after correction, which is an indication of the poor performance of the GACOS approach. In this example, the phase elevation method performs better as the signal is concentrated on the high relief areas (Figure 3c), and the standard deviation decreases to 0.9 cm after correction (Figure 3f). In the same way as the previous example, neither GACOS or phase elevation approaches fit the data at high elevations (above 2,200 m) (Figure 3d), which results in a small positive residual on Agung's summit.

Selected Examples
From these two examples, it is evident that a single solution may not be suitable for all cases and that even in a single location, different approaches are suited to the atmospheric conditions on different days. In particular, the GACOS corrections perform well for long wavelength and large-magnitude signals whereas the phase elevation could be more appropriate for subtle tropospheric signals around volcanic edifices. In the next section, we consider the entire data set and evaluate the overall quality of these two atmospheric corrections using a statistical approach.  propose four metrics for assessing the performance of GACOS corrections. Two of these depend solely on location and are constant through time: (1) the time difference between satellite acquisition and ECMWF model and (2) the local relief. These metrics are useful for comparing the potential benefits and challenges of using GACOS corrections at different volcanoes. In the case of Agung, the ascending track was acquired at 5:41 p.m. and the descending track at 4:53 a.m. Jakarta local time. The corresponding time difference to the HRES-ECMWF model outputs is 1 hr 19 min and 2 hr 7 min, respectively, compared to a maximum value of 3 hr. Given the rapidly varying properties and high relief associated with small tropical islands, these values suggest that Agung will be a challenging target for GACOS, and hence is a good test for studies of global applicability.

Performance Metrics
The other two metrics depend on the atmospheric conditions at the time of acquisition: (1) the reduction in the standard deviation (std) after the correction and (2) the internal cross validation of the GACOS correction using the ECMWF grid. These could be used for flagging unusual atmospheric conditions within a time series of data from a single track and are described further in this section.

Reduction in Variability
For each interferogram, we estimate the reduction in standard deviation defined by Δ = i − i , where i is the standard deviation of the uncorrected interferogram and is the standard deviation of the corrected one (Tables S1 and S2 in the supporting information). Uncorrected interferograms generally have a higher standard deviation on the ascending track (4.5 ± 1.8 cm) than the descending track (4.2 ± 1.5 cm), which likely relates to acquisition time. In tropical regions, we expect more turbulent atmospheric signals in the afternoon due to solar heating than in the early morning before sunrise. For Agung, the ascending track is acquired in the early evening (5:41 p.m., Jakarta local time), whereas the descending is acquired in the early morning (4:53 a.m, Jakarta local time).
We assume that atmospheric signals dominate the short-period (12 days) interferograms such that positive values of Δ indicate a reduction in atmospheric artifacts, while negative values indicate that the correction added noise to the data. Figure 4 shows the parameter Δ (in %) plotted as a function of the uncorrected standard deviation i for the two approaches.
For the GACOS method, the corrections are very efficient for large atmospheric signals ( i > 2 cm) and Δ can reach 25% to 50% (strong signal in Figure 4a), as shown by one example in Figure 2. However, if the original variability is small ( i < 1.5 cm), corresponding to calm atmospheric conditions, GACOS model produces an inaccurate correction with negative values of Δ (low signal in Figure 4a), as illustrated by the example shown in Figure 3. The trend shows that the quality of GACOS corrections is highly sensitive to the amplitude of the atmospheric signal. Assuming that the precision of the weather models remains constant over the data set, the GACOS prediction is more likely to be very different from data values for small tropospheric signals. When the amplitude of the atmospheric signal is larger than the model precision, the discrepancies between data and model are reduced and by consequence the quality of corrections improves.
For the phase elevation method, the reduction in standard deviation, Δ , is generally low: only 26% of ascending interferograms and 13% of descending interferograms have Δ values above 10% (Figure 4b). Compared to GACOS corrections, Δ values cannot be negative as the empirical method only removed correlated signals; however, we cannot be certain that the signal removed is only atmospheric and not from other sources. There is no significant trend (the coefficient of determination, R 2 ∼ 0.2) and Δ is larger than 25% only for few ascending dates. Using this metric alone, the phase elevation correction appears to perform poorly; however, visual inspection of individual examples suggests that it can be advantageous in certain circumstances (e.g., Figure 3). Overall, the GACOS corrections perform better than the phase elevation model, especially for strong atmospheric signals associated with large original standard deviations. However, even the HRES-ECMWF products do not have sufficient spatial resolution to capture small-scale atmospheric signals (<10 km).
We therefore consider strategies combining the benefit of the two methods with the aim of mitigating a wider range of atmospheric artifacts. By only applying the GACOS corrections, we notice that in a few interferograms the corrected phase still shows a strong correlation with the topography, such as the interferogram 21 May 2017 to 2 June 2017 ( Figure 5). In this example, the phase is highly correlated with topography (Figures 5a and 5b), but the GACOS correction does not capture the pattern and it only removes a small portion of the atmospheric signals located along the northern coastline ( Figure 5d). The reason could be that the spatiotemporal variability of the atmosphere was too high at the time of the acquisition to be captured by the ECMWF 6-hr model. After GACOS corrections, the correlation with topography is even clearer with a large positive signal along the entire volcanic chain (Figure 5f). By applying the phase elevation correction after the GACOS corrections (Figure 5e), it is possible to remove the correlated residual signal and obtain an interferogram with almost no residual signals (Figure 5g). The coefficient of correlation between the data and the correction, R data , reaches about 0.8 for the combined approach, which is twice the value obtained for GACOS correction. For comparison, the reduction of the standard deviation is 2.9, 23.5, and 32.3% for the GACOS, the phase elevation and the combined corrections, respectively. This case illustrates how the combined approach can be significantly more efficient than either individual correction.
Looking at the overall data set, the phase elevation approach does not produce a significant reduction in the global standard deviation ( Figure 6) as it only removes stratified signals correlated with the relief of Agung and Batur, while other large-scale atmospheric features remain uncorrected (Figures S1b and S2b in the supporting information). In contrast, the GACOS corrections produce a large reduction in standard deviation for most of the interferograms (Figures S1c and S2c in the supporting information), and particularly for peak values in July 2017 or March 2018 for ascending data (Figure 6a) and May or June 2017 for descending data (Figure 6b). The combined approach discussed above produces further corrections for some interferograms and provides the lowest standard deviation values for most of the interferograms (Figure 6, blue line).
Although the reduction of standard deviation with the combined approach is similar to the GACOS values for our data set, the benefit is evident for individual interferograms ( Figure 5). In the framework of an early-warning system, it is crucial to obtain interferograms almost free from atmospheric signals to support a fast response and reliable interpretation. We believe the implementation of combined corrections is and (e,g) correction and residual using a combined approach (1 = GACOS correction followed by 2 = a phase elevation correction, see subset). The coefficient of correlation between the correction and the data, R data , is indicated on the bottom left corner. The standard deviation (in radian) of each corrected interferogram, , is shown on the bottom right corner. appropriate as it reduces the possibility of misinterpretation of tropospheric signals and ground deformation signals, such as that which occurred during the 2017 unrest at Agung.

Internal Validation
In addition, we use the cross-validation approach to assess the interpolation performance of ZTD maps produced by GACOS. Cross validation has the advantage of being available a priori, before the interferogram is even processed. Cross validation is conducted as follows: one grid point from the regular 0.125 • spacing grids was excluded and the ZTD values from the other grid points were used to determine the ZTD value at the grid point considered. This procedure was repeated for each grid point and the root-mean-square (rms) difference computed between the interpolated and original ZTD values. Because of the small and regular spacing (and therefore high spatial correlation) of the HRES-ECMWF ZTD data points, the cross validation was conducted using 20% of the grid points chosen at random to reduce the spatial correlation. In theory, the most precise GACOS delays maps will be associated with the lowest rms values.
For each interferogram, Figure 7 shows a comparison between the rms values from cross validation (black lines) and the standard deviation of the corrected interferograms ( ) for ascending (blue line) and descending tracks (red line). All time series have been normalized using X−X min X max −X min to range between 0 and 1 and make the comparison between the two performance metrics easier. For the ascending data, we observe that some high rms values are correlated with high std values (Figure 7a) (e.g., point #11, interferogram 6 September 2017 to 18 September 2017). Such cases occur when both metrics are consistent and it confirms that the performance of the GACOS correction will be poor for these dates. However, the overall correlation between the two metrics remains low with R 2 = 0.09 (Figure 7a). For the descending data, the correlation is even weaker with R 2 < 10 −5 , which suggests that the two metrics are independent (Figure 7b). These examples highlight the limitations of the cross-validation metric: in some cases, it predicts a good quality GACOS correction, whereas the corrected interferogram still contains turbulent signals causing a high standard deviation.
In previous studies focusing on California and United Kingdom, the cross-validation metric showed the highest potential for assessing the performance of the GACOS atmospheric corrections . However, in these cases, the cross-validation metric was calculated using GPS measurements whereas at Agung, there is no GPS network and the cross-validation scores are only calculated using the ECMWF grid,  which is already spatially smoothed. Without a dense GPS network, the cross-validation metric does not capture the poor performance of the delay maps when small wavelength turbulent signals are common. Thus, we conclude that the cross-validation metric is less well suited for tropical island volcanoes than the reduction in variability. As a consequence, it is not possible to predict when corrections will work well a priori, and it is necessary to apply the GACOS correction to the data before its performance can be accurately assessed.

Detecting Deformation Anomalies
The use of satellite data was critical for monitoring ground deformation during the 2017 seismic swarm at Agung. Due to the presence of strong atmospheric signals associated with steep tropical islands, the interferograms were initially misinterpreted. Even once the atmospheric signals had been identified, the threshold for detecting deformation was very high and confidence was low. Consequently, the characterization of the preeruptive deformation and the modeling of the source(s) was performed after the eruption . To be useful in volcanic crises, monitoring systems need to be available quickly, robust, and objective. Automated analysis tools are particularly useful for volcano observatories as they can be routinely applied at large number of systems and used to flag anomalous events for expert analysis (Anantrasirichai et al., 2018(Anantrasirichai et al., , 2019Gaddes et al., 2018).

Routine Production of InSAR Time Series
We hypothesize that near-real-time production of corrected InSAR time series would have reduced uncertainty regarding the cause of the seismic swarm, improving the quality of scientific advice provided to the Indonesian authorities who were managing the crisis. In this section, we first compare the displacement time series from May 2017 to May 2018 with and without atmospheric corrections (Figures 8 and 9), and then assess automatic algorithms for detecting deformation based on the analysis of the Agung time series (Figure 10). Figure 8 shows the cumulative LOS displacements at the time period of the seismic swarm for the ascending and the descending tracks, respectively. For the cases without atmospheric correction (Figures 8a and 8c), the phase difference is still dominated by large wavelength signals due to atmospheric delays, which are partially masking the ground deformation signal. For the cases with combined atmospheric corrections (GACOS + phase elevation), an uplift signal is detected in early October on the northern flank of Agung for both ascending and descending tracks (Figures 8b and 8d). The difference between the spatial extent of the signal between ascending and descending tracks is due to the difference in LOS direction . Figure 9 shows the time series of cumulative LOS displacements at two locations: point A, which is the center of the uplift signal detected in the cumulative ground deformation maps and point B, which is the center of Batur caldera, considered to be a nondeforming reference location (Figure 8). The uncorrected time series (Figures 9a and 9c) show large variations (±5 cm) due to atmospheric artifacts. No clear ground deformation trend can be extracted from the ascending track, and the deformation signal on the descending track would not have been convincing alone. In the corrected time series (Figures 9b and 9d), the amplitude of the atmospheric noise is significantly reduced. Both ascending and descending track data show a clear step change at the time of the seismic swarm, 1 month before the eruption.

Automatic Detection of Unrest Period
Here, we test two methods for the automatic detection of rapid unrest: (1) the simple thresholding of the time series and (2) a sequential anomaly method. The performance of both methods will be evaluated using the Receiver Operating Characteristic (ROC) curves. ROC curves test the ability of a classifier system to discriminate between two classes (here, unrest and no unrest) by plotting the true positive rate (e.g., probability of detection) against the false positive rate (e.g., probability of false alarm) at various threshold. The Area Under the Curve (AUC) quantifies the success of the classification with 0.5 being a worthless test (the system has no ability to distinguish between the two classes) and 1 being a perfect test (100% detection with no false alarm). Based on the LOS time series (Figure 9), we consider that the unrest period lasted between October and mid-November 2017 .

Simple Threshold
One approach to detecting rapid changes is to apply a simple threshold to the incremental time series (Figure 10). However, the detection performance is controlled by the balance between the quality of the data and the threshold value, which needs to be selected for each system. Figures 10a and 10b shows the uncorrected and corrected incremental time series at the center of uplift (point A in Figure 8). Without atmospheric corrections, we do not observe significant anomalies during the unrest period (Figure 10a). With atmospheric corrections, we start to distinguish some high values during the unrest period (Figure 10b). The limitation of the method is that we have no a priori information about the optimal threshold value to chose, as it depends on too many factors (e.g., quality of the time series, atmospheric conditions, or amplitude of the ground deformation). As a result, a low threshold will increase the number of false alarms (false positive) whereas a high threshold will increase the number of missed detections (false negative).
To evaluate the performance of the detection method, we build the ROC curves by calculating the True Positive Rate ( ascending and descending time series (Figure 11a), and values below 0.5 confirm that the positive and negative classes cannot be differentiated. After atmospheric correction, the performance of the thresholding algorithm improves, but the AUC values are still 0.78 and 0.82 for ascending and descending time series, respectively, indicating a ∼80% chance the approach can distinguish between positive (unrest) and negative (no unrest) classes (Figure 11b). A high TPR is only achieved when the detection threshold is small (0.5-1 cm) leading to numerous false alerts (FPR of 0.7 and 0.4 for ascending and descending time series, respectively).

CUSUM Control Chart
Given the challenges associated with detecting anomalies using simple thresholds in InSAR time series, we test a method for detecting consecutive anomalies based on the fact that deformation tends to be consistent with time, whereas anomalies caused by unusually strong atmospheric artifacts can be identified by a sign reversal in sequential interferograms (Massonnet & Feigl, 1995). We adopt the cumulative sum (CUSUM) control chart, which is a sequential analysis technique developed by E. S. Page of the University of Cambridge (Page, 1954), typically used for monitoring small changes in the process mean. In Earth Sciences, the method has been successfully applied to detect land cover changes from MODIS time series (Grobler et al., 2012;Kucera et al., 2007;Rahman et al., 2013).
The CUSUM method tests if the new value of a time series is statistically different from the previous samples using the cumulative sum of deviations from the target mean. Two criteria are calculated, U i and L i , which track positive and negative sequential anomalies, respectively. 10.1029/2019JB017908 where z i is the z-score of the sample i (i.e., the number of standard deviations from the mean value) and n is the number of standard deviation required for the shift to be detectable. The choice of n = 1 corresponds to the detection of one-sigma shifts in the mean. For each pixel, we calculate the sum of U i +L i at time t i based on the mean and standard deviation of all the previous samples (t < t i ). This simulates the situation where the detection is performed in real time only with the past temporal information. The method is cumulative as the contribution is summed over time such that a sequence of values statistically different from the background noise will have higher U i (or L i ) values than a single-point anomaly. Figure 10c shows the CUSUM results applied on the corrected incremental time series, using n = 1. For the descending time series, U i + L i = 2.2 on 3 October and U i + L i = 2.8 on 15 October (Figure 10c, red), corresponding to confidence levels of 97% and 99.5%, respectively that the data on these dates are statically different from previous samples. The ascending track is noisier than for the descending track due to poorer atmospheric corrections and therefore the anomalies are smaller in amplitude with U i + L i = 1.5 on 12 October and U i + L i = 1.9 on 5 November (Figure 10c, blue), equivalent to a probability of 87% and 94%.
Other single peaks can be identified through the time series such as in July 2017 (U i + L i = 1.4) for the descending track or in September 2017 (U i + L i = −0.8) for the ascending track.
In the CUSUM control chart, deformation signals will appear as high values on consecutive dates, under the conditions that the duration of the unrest is larger than the temporal sampling (12 days here). In contrast, atmospheric residual signals are not correlated in time and show as a single peak that reverses in sign on the next date. The joint evaluation of ascending and descending CUSUM values on consecutive dates improves the flagging of the period of unrest. Based on this criterion, the CUSUM flowchart shows highest values for the period between October and mid-November (Figure 10c). The sum of ascending and descending anomalies exceeds 2, which indicates a level of confidence above 95%. This time period corresponds to the ground deformation associated with the dike intrusion between Mount Agung and Batur identified retrospectively by Albino et al. (2019). The other anomalies detected in the time series are single events (e.g., present in a single date and/or only in one track), which correspond to residual atmospheric noise (Figure 10c).
We calculate the ROC curves by calculating the FPR and the TPR at different threshold values between 0 and 3 . The CUSUM approach outperformed the threshold method with AUC values of 0.95 and 1 for ascending and descending, respectively ( Figure 11c). Because the parameter n (equation (1)) could affect the performance of the detection, we calculate the AUC values for different values of n ( Figure 11d): (i) using high values (n > 1) decreases the performance by increasing the number of missed detections (e.g., only large amplitude changes are detected); (ii) using small values (n < 0.5) decreases the detection performance by increasing the number of false detections (e.g., more small changes associated with noise are detected). Our initial choice of n = 1 is justified as it falls in the range of values providing the maximum performance (gray area in Figure 11d).

Toward Real-Time InSAR Monitoring of Volcanic Ground Deformation
The 2017 volcanic crisis at Mount Agung showed once again the potential of InSAR to identify ground deformation signals and constrain the geometry and depth of the plumbing system . However, the robust detection of the preeruptive uplift was delayed due to the presence of strong and variable atmospheric artifacts in the data set and therefore the origin of the seismic swarm was revealed only after the beginning of the eruption. This example points out the importance of real-time InSAR analysis in providing ground deformation monitoring.
Here, we propose a new processing flow chart for the near-real-time automated detection of volcanic deformation ( Figure 12). New interferograms are calculated and the time series are updated each time a new Sentinel-1 SAR image is acquired over the area of interest. The innovative aspect of our processing consists of two additional toolboxes (shaded boxes in Figure 12), which contribute to improving the identification of InSAR ground deformation signals: (a) the systematic corrections of tropospheric delays based on combination of high-resolution ECMWF weather models and phase elevation correlation and (b) the statistical analysis of corrected time series to detect sequential anomalies. Our approach allows the fast production and analysis of InSAR time series (within 2 days), which will support volcano observatories (e.g., the Indonesian Center for Volcanology and Geological Hazard Mitigation) for eruption forecasting and risk mitigation. Processing flow chart proposed to automatically detect in near-real-time volcanic unrest using Sentinel-1 SAR data. New aspects consists of (1) the systematic correction of tropospheric delays combining high-resolution weather models and phase elevation corrections and (2) the detection of sequential anomalies in the InSAR corrected time series using cumulative sum control chart. The duration of each step is indicated by bold numbers.
Because our methodology only requires open-access data available online (see dark rectangles in Figure 12), it can be applied to other volcanic regions. As stated in section 1, there are several challenges for detecting subtle deformation signals on andesitic stratovolcanoes using InSAR (Pinel et al., 2011;Parks et al., 2011;Remy et al., 2015). Our automatic algorithm partially addresses these limitations. This approach is suitable for the detection of small rapid changes in stationary time series and performs better for cases with a high signal over noise ratio. However, it is less well suited to the analysis of long-duration trend signals observed at large silicic systems such as Campi Flegrei (Italy), Corbetti (Ethiopia), or Uturuncu (Bolivia) (Henderson & Pritchard, 2013;Lloyd et al., 2018;Trasatti et al., 2008).

Conclusions
Our study shows the importance of correcting atmospheric delays on InSAR data at tropical volcanoes, using the example of 2017 volcanic crisis at Mount Agung (Bali, Indonesia). The quality assessment of atmospheric corrections over 62 interferograms covering 1-year period (May 2017 to May 2018) demonstrates that high-resolution weather models have a better performance than the simple phase elevation approach. However, the spatial resolution of weather data remains too coarse to accurately correct some of the small atmospheric signals located at the top of high peaks.
After correction, the InSAR time series show a wide preeruptive uplift between Agung and Batur volcanoes during October 2017, which correlates in time with the peak of the seismicity. Our anomaly detection approach enables us to detect unrest at Agung with 95% confidence only 2 weeks after the start of the seismic swarm and 1 month before the eruption. Despite its simplicity, our approach can flag anomalous behavior in near real-time (within 2 days of acquisitions), without requiring expert interpretation, and hence shows potential as a monitoring tool for volcano observatories globally.