The Application of Convolutional Neural Networks to Detect Slow, Sustained Deformation in InSAR Time Series

Automated systems for detecting deformation in satellite interferometric synthetic aperture radar (InSAR) imagery could be used to develop a global monitoring system for volcanic and urban environments. Here, we explore the limits of a convolutional neural networks for detecting slow, sustained deformations in wrapped interferograms. Using synthetic data, we estimate a detection threshold of 3.9 cm for deformation signals alone and 6.3 cm when atmospheric artifacts are considered. Overwrapping reduces this to 1.8 and 5.2 cm, respectively, as more fringes are generated without altering signal to noise ratio. We test the approach on time series of cumulative deformation from Campi Flegrei and Dallol, where overwrapping improves classification performance by up to 15%. We propose a mean‐filtering method for combining results of different wrap parameters to flag deformation. At Campi Flegrei, deformation of 8.5 cm/year was detected after 60 days and at Dallol, deformation of 3.5 cm/year was detected after 310 days. This corresponds to cumulative displacements of 3 and 4 cm consistent with estimates based on synthetic data.


Introduction
Interferometric synthetic aperture radar (InSAR) can be used to measure ground displacement over large geographic areas.However, while detecting deformation using InSAR images is conceptually straightforward, it is difficult to automate.For volcanic regions, atmospheric water vapor is a particular challenge since the stratification of atmospheric water vapor forms a topographically correlated pattern that can be difficult to distinguish from the deformation caused by the pressurisation of magma reservoirs beneath volcanic edifices (Beauducel et al., 2000;Ebmeier et al., 2013;Remy et al., 2015).In urban environments, the deformation patterns typically have a lower rate and smaller spatial scale but are more distinct from atmospheric noise.
Previous studies have shown that deep convolutional neural networks (CNNs) have the capability to identify volcanic deformation signals in wrapped interferograms when trained on large data set (Anantrasirichai et al., 2018(Anantrasirichai et al., , 2019;;Valade et al., 2019).The high-frequency content of wrapped fringes provides strong features for machine learning algorithms, and Anantrasirichai et al. (2018Anantrasirichai et al. ( , 2019) ) use a method based on edge detection.The outputs are expressed as a probability, which can be used to flag deformation.However, for C-band satellites, one fringe corresponds to 2.8 cm of deformation or 85 cm/year for a 12-day interferogram and 170 cm/year for 6-day interferogram.The approach of (Valade et al., 2019) also only tested short-term interferograms that showed deformation of >10 cm Global Volcanism Program (2019).These high rates are typically only observed for very short periods associated with dike intrusions or eruptions (Biggs & Pritchard, 2017).Yet there are many deformation signals that occur at lower rates but for longer duration, such as sustained uplift at silicic volcanoes (Henderson & Pritchard, 2017;Lloyd et al., 2018;Montgomery-Brown et al., 2015;Remy et al., 2014;Trasatti et al., 2008), subsidence and heave in former coalfields (McCay et al., 2018), engineering projects such as tunnelling, and landsliding of natural and engineered slopes (Whiteley et al., 2019).
Here, we investigate two adaptations to enable a machine learning system to detect slow, sustained deformation: 1) the use of cumulative time series of displacements to increase the time interval and hence the signal to noise ratio and 2) rewrapping at different gain each time displacement map to generate additional

CNNs
CNNs are deep feed-forward artificial neural networks (Goodfellow et al., 2016), designed to take advantage of 2-D structures, such as an image.In this study, we use a transfer-learning strategy, employing the fine-tune AlexNet (Krizhevsky et al., 2012) developed by Anantrasirichai et al. (2019).The wrapped interferograms are employed because the high-frequency content at fringes is easy to identify and provides strong features for distinguishing deformation from other signals.The wrapped interferogram is first converted into a grayscale image (i.e., the pixel values are scaled to [0,255]) and divided into overlapping patches at the required input size for AlexNet (224 ×224 pixels).Following Anantrasirichai et al. (2019), the top-left position of each patch is then repeatedly shifted by 28 pixels to cover the entire image.The output of the prediction process is a probability P of there being deformation in each patch, and the probabilities from overlapping patches are merged using a rotationally symmetric Gaussian low pass filter with a size of 20 pixels and standard deviation of 5 pixels.The network was initially trained with synthetic InSAR data representing deformation, turbulent, and stratified atmospheric contributions and then retrained with a combined data set consisting of synthetic models and selected real examples (the accuracies of the networks trained by the synthetic data set alone and the combined one are 98.94% and 99.96%, respectively Anantrasirichai et al., 2019).The interferogram is flagged as containing deformation when P > 0.5.

InSAR Dataset
Initially, we test the detection threshold of the network using synthetic data.We create a data set of 405 synthetic interferograms (X) using three components namely deformation D, stratified atmosphere S, and turbulent atmosphere T, using the linear function X = D + S + T, where ,  ∈ {0, 0.5, 1}.Synthetic deformation signals, D, are generated using a point pressure source (Mogi) model (Mogi, 1958), which demonstrates surface deformation associated with inflation and deflation of a magma chamber.In this study, D is modelled with depths of 3, 4, and 5 km, incidence angles of 1 • , 23 • , and 44 • , and volumes of 10  m 3 , where  ∈ [5, 5.5, 6, … , 7].For stratified atmosphere S, we use the signal at Beru volcano (4 January to 17 March 2017) via the Generic Atmospheric Correction Online Service.It is based on weather model data (Yu et al., 2018), and a zenith total delay map is derived from the high-resolution water vapor delays generated by the European Centre for Medium-Range Weather Forecasts.Turbulent atmospheric delays, T, are spatially correlated, and their covariance can be described using an exponentially decaying function.The one-dimensional covariance function is c i =  2 max e (−d i ) , where c ij is the covariance between pixels i and j, d ij is the distance between the pixels,  2 max is maximum covariance, and  is the decay constant, which is equivalent to the inverse of the e-folding wavelength (Biggs et al., 2007).We employ  2 max of 7.5 mm 2 and  of 8 km (Anantrasirichai et al., 2019).
We then test our approach on two examples: Campi Flegrei (Italy) where a new pulse of uplift began in July 2017 and Dallol (Ethiopia), which is uplifting at a steady rate of 3.5 cm/year (Figure 3).Campi Flegrei is an active 13-km wide caldera located within the bay of Naples, a densely populated area of 3 million inhabitants (Global Volcanism Program, 2013).Since the last eruption in 1538 at Monte Nuovo, the caldera system showed signs of unrest with episodes of rapid ground uplift (e.g., 1950-1952, 1969-1972, and 1982-1984) followed by long period of recovery associated with slow ground subsidence (e.g. 1984-2010;Del Gaudio 2010;Troiano, 2011.Dallol is an Ethiopian volcano located in the Danakil depression.The crater lies 48 m below the sea level and was formed during a phreatic eruption in 1926 (Global Volcanism Program, 2013).Recent unrest associated with the emplacement of a magma intrusion have been detected by InSAR in 2004 (Nobile et al., 2012), and a phreatic eruption have been recorded on January 2011 (Global Volcanism Program, 2013).
For each volcano, we process Sentinel-1 SAR images using the LiCSAR system (González et al., 2016), which generates automatically the three short-duration interferograms for each time acquisition (i.e., 6, 12, and 18 days for Campi Flegrei and 12, 24, and 36 days for Dallol) with a spatial resolution of ∼100 m (Li et al., 2016).Interferograms are not corrected from atmospheric signals as our CNN network was initially trained to distinguish between ground deformation and atmospheric signals.A linear least-square inversion is performed on the network of unwrapped interferograms to obtain the time series of cumulative ground deformation using singular value decomposition and assuming no deformation at the first date (Berardino et al., 2002;Schmidt & Bürgmann, 2003;Usai, 2003).

10.1029/2019GL084993
Figure 1.Example synthetic interferograms rewrapped with  = 1, 2, 4, and 8 (left column) and with shifted wrap boundaries with  = 0, /2, , and 3/2 (right column).The relative weighting of the stratified and turbulent atmospheric components are given by  and , respectively.Top row = = 0, middle row = = 0.5, and bottom row  = = 1.The probability, P, that the image contains a deformation signal according to the CNN is listed beneath each image.

Data Wrapping
The phase  of an interferogram varies between − and  and is typically unwrapped before geophysical analysis (Chen & Zebker, 2002).CNNs have been previously applied to wrapped InSAR data (Anantrasirichai et al., 2018;Valade et al., 2019), which is preferable for detecting rapid deformation signals as the unwrapping process is time consuming and can introduce errors.However, for studying slow, sustained deformation, unwrapping is necessary to produce time series of cumulative deformation.Nonetheless, fringes provide strong features for machine learning and in this study, we rewrap each displacements map deduced from the least-square inversion and use the pretrained network of Anantrasirichai et al. (2019).
We test whether altering the number and location of the wrap boundaries has the potential to improve the ability of the CNN to detect deformation.Reducing the wrapping interval increases the number of fringes without altering the signal to noise ratio.We generate a new phase  ′  with a wrap gain  using ' ≡  (mod 2), where  is a positive integer.The wrap interval is reduced by 1∕ of the original phase value, and the number of fringes increases  times, for example, = 2 produces twice as many fringes (Figure 1; left).We also consider shifting the wrapping boundaries by adding a constant phase offset  to the unwrapped interferogram prior to rewrapping using  ′  , that is, ' ≡ + (mod2).The phase discontinuities occur in physically arbitrary locations, so for some cases, the number of fringes will increase but in others it will decrease (Figure 1; right).

Detection Thresholds
Figure 1 (left) shows examples of synthetic interferograms with wrap gains = 1, 2, 4, and 8.For each example, the CNN output probability is below the detection threshold (P < 0.5) at  = 1 but greater than the detection threshold (P > 0.5) at = 2, 4, and 8.This illustrates that the increased number of fringes improve the ability of the CNN to detect the deformation, even though the signal to noise ratio of the inteferograms remains constant.In contrast, shifting wrap boundaries does not always improve the detection performance as shown in Figure 1 (right), where P may increase and also decrease.
Figure 2 shows the output probability, P, against maximum displacement for each interferogram in the synthetic data set.The transition between P ∼ 0 (undetectable) and P ∼ 1 (detectable) occurs over a narrow 10.1029/2019GL084993 Figure 2. Analysis of the detection threshold of the convolutional neural network using synthetic data with a range of source models and atmospheric conditions.The probability, P, that the image contains a deformation signal according to the convolutional neural network is plotted against the maximum displacement for a range of simulations.The range of deformation signals, D, is created using a Monte Carlo approach and a point (Mogi) model.The deformation is combined with stratified (S) and turbulent (T) atmospheric components according to the equation X = D + S + T. and  are set equally to 0, 0.5, and 1 in the first, second, and third rows, respectively, and more combinations of  and  are shown in the supporting information (Figure S2).The last row shows the detection thresholds when P= 0.5 for different values of  and .
range and corresponds to the minimum value of displacement that the CNN can identify.We fit a sigmoidal curve, defined as f(x)=(1+e −a(x−b) ) −1 , to estimate the deformation detection threshold for each of the different levels of stratified noise () and turbulent noise ().As expected, stronger stratified atmosphere (increasing ) and stronger turbulent atmosphere (increasing ) increase the detection threshold from 3.9 cm for  =  = 0 to 6.3 cm for  =  = 1.Conversely, increasing the wrap gain from  = 1 to  = 2 reduces the detection threshold by 25-30%, for example, for  =  = 0.5 and 1, giving final values of 3.8 and 5.2 cm, respectively.However, although the higher wrap gains produce a lower detection threshold overall, there is significant scatter in the probabilities, suggesting that such a high wrap gain would produce a large number of false-positive results.Shifting the wrap boundaries produces a small reduction in detection threshold for noise-free data but no change in the overall detection threshold when noise is included in the simulation as shown in Figure 2 (bottom row).
The estimated detection thresholds can be applied either to a single interferogram or to a cumulative time series, where they correspond to the deformation rate multiplied by the time-series duration.For example, the detection threshold of 5.2 cm for a volcanic environment corresponds to a steady rate of 1.3 cm/year for the 4 years of Sentinel-1 data available at the time of writing (2015)(2016)(2017)(2018)(2019).This is a significant improvement 10.1029/2019GL084993 to the minimum rate of just 1.80 m/year reported in Anantrasirichai et al. (2018) for applying CNNs to individual 12-day wrapped interferograms.

Application to Real Examples
For both the Campi Flegrei and Dallol data sets (Figure 3), the CNN outputs a probability map for each step of the time series, and we plot the probability values at the stable point A (Figure 3; third row) and the volcanic center, B (Figure 3; fourth row).For each example, we test four wrap gains, that is, = 1, 2, 4, and 8.
For Campi Flegrei (the deformation rate of 8.5 cm/year), the probability at point A is <0.05 for almost all time steps, whereas at point B, there is an increase in P around the start of the deformation.Using the original wrap gain ( = 1), we detect deformation in February 2018, 7 months after the episode began.At this point, the cumulative displacement is > 5 cm, consistent with the detection threshold in our synthetic tests.In the times series with  equal to 2, 4, and 8, we identify deformation at ∼3 cm (2.5 months from onset), ∼2 cm (1.5 months from onset), and <1 cm (−9 months from onset), respectively.However, the higher wrap gains ( > 4) incorrectly exceed the detection threshold of P = 0.5 before the deformation began, with a total of 15 false positives for  = 8.At Dallol, the deformation time series is noisier, and the rate of deformation (3.5 cm/year) is lower than the previous example.At  = 1, the probability first exceeds P = 0.5 after 21 months but continues to drop below the threshold throughout the entire 4-year time series (Figure 3; right).The cumulative deformation after 4 years is 10 cm, significantly greater than the synthetic detection threshold.Visual inspection shows that the low probabilities occur when the surrounding area is incoherent, a situation that was not included in the synthetic tests.In these cases, doubling the wrap gain to  = 2 increases the probability above the 10.1029/2019GL084993 threshold (P = 0.5).The highest wrap gain =8 identifies deformation at point B in only the second time step when the displacement is < 1 cm, but it falsely identifies deformation at point A for the first 9 months.Examples of wrapped interferograms at Campi Flegrei and Dallol when  = 1 and 2 are shown in Supporting Information Figure S3, and the interferograms causing false positives at Dallol are shown in Supporting Information Figure S4.
Increasing the wrap gain () of the interferograms reduces the detection threshold for the CNN thus allowing us to detect slow, sustained deformation earlier in time series.However, increasing the wrap gain too much causes the CNN to misidentify features caused by atmospheric artifacts.Consequently, decisions based on only large values of  would have many false positives while those using only small values of  might not be able to detect slow deformation.We therefore compute the final probability P by combining the probability at a range of wrap gains (equation ( 1)), using N = 4. (1) The value of P shows a similar trend to the magnitude of displacement -low P at small displacement and large P at large displacement (Figure 3; last row).The probability first exceeds the threshold P = 0.5 on 27 August 2017 for Campi Flegrei and 2 December 2015 for Dallol (Figure 3), corresponding to 2 months after the onset of deformation and 11 months after the start of the time series, respectively.

Discussion
We evaluated the results from the examples at Campi Flegrei and Dallol using a receiver operating characteristic (ROC) curve (Figure 4), where true positive and false positive rates were calculated by varying the probability thresholds to identify deformation and non-deformation.The true positive rate is the fraction of predicted deformation that are retrieved over the total number of actual deformation, while the false positive rate is the number of nondeformation wrongly identified as deformation divided by the total number of actual nondeformation.We computed the area under the ROC curve (AUC): good classifiers will give high AUC values as they detect the positive signals correctly and few true negatives are falsely identified.For Campi Flegrei, the highest AUC (0.989) is produced by  = 6, which indicates good separation between classes.The AUC decreases for higher wrap gains as the number of false positives increase.For Dallol, the ROC curves are computed assuming the deformation starts on the first date of the time series (15 January 2015), and the highest AUC results (0.985) are for wrap gains of = 7. The AUC under the curve for the combined probability P is 0.989 and 0.984 for Campi Flegrei and Dallol data sets, respectively.This demonstrates that CNNs can possibly be applied to wrapped InSAR data to detect deformation at rates as low as 3.5 cm/year.
Here, we have tested the ability of CNNs to detect slow, steady deformation in volcanic environments, but further testing is still required to assess whether this technique could also be used to detect anthropogenic sources of deformation, such as subsidence associated with resource extraction or infrastructure instability.The characteristics of the deformation are quite different but so are the noise characteristics.In particular, the deformation has a smaller spatial extent and the topographic gradients smaller, which means that it is possible to reduce the atmospheric components by filtering (Ferretti et al., 2011;Hooper et al., 2004).Thus, although the rates of deformation are lower, we might expect similar signal-to-noise ratios.
Our approach is based on the ability of CNNs to distinguish between different sources using the pattern of fringes (Anantrasirichai et al., 2018(Anantrasirichai et al., , 2019)).However, to create cumulative time series, the interferograms need to be unwrapped first, and we artificially rewrap the data.Alternatively, CNN-based methods applied directly to unwrapped data may also have the potential to detect slow, steady deformation, and a direct comparison of the performance of the two approaches using robust performance metrics (e.g.precision and recall) on a suitable range of case studies would be interesting.

Conclusions
In this paper, we study the feasibility of using CNNs to detect slow surface deformation in InSAR images by rewrapping cumulative time series.We first tested the technique on synthetic interferograms based on volcanic environments and found that when the number of fringes are doubled, the detection threshold is lowered by 25-30%.The detection threshold of 5.2 cm corresponds to rates of 1.3 cm/year for the 4-year 10.1029/2019GL084993 Sentinel-1 time series currently available.Subsequently, we test using the cumulative times series from two volcanoes: Campi Flegrei and Dallol and demonstrate the ability to detect rates as low as 3.5 cm/year, with an improvement of the AUC up to 15% when overwrapping is applied.However, increasing the wrap gain too high increases the number of false negatives, so we propose a method to combine the results from multiple wrapping gains.This shows potential for future use, though further adaptations may be required to apply the technique to a wider range of environments.

Figure 3 .
Figure 3. Examples using Sentinel-1 time series from Campi Flegrei (left column) and Dallol (right column).Top-row figures are the unwrapped line-of-sight cumulative displacements at the end of period of study with the locations of two studied areas, A and B, and the reference point (REF).Second-row figures show the time series of cumulative displacements.The ground at area A is considered stable, whereas B is located at the volcanic center showing significant deformation.The CNN-output probability for points A and B are shown in the third and fourth rows for wrap gains = 1,2,4, and 8.The fifth row contains the average probability from the selected wrap gains as defined by equation (1).The green and the red dashed lines correspond to the deformation detection threshold (P = 0.5) and the start date of deformation, respectively.

Figure 4 .
Figure 4. Performance analysis of the convolutional neural network approach.The top row shows the receiver operating characteristic curves calculated from both points Aand B, while the bottom row demonstrates the area under the receiver operating characteristic curves for wrap gains  = 1,2,4, and 8.The higher area under the receiver operating characteristic curve indicates better performance (fewer false positives and false negatives).(Left column) Campi Flegrei and (right column) Dallol.The start date of deformation is taken to be 1 July 2017 for Campi Flegrei and 15 January 2015 for Dallol, which corresponds to the start of the time series.