A New InSAR Phase Demodulation Technique Developed for a Typical Example of a Complex , Multi-Lobed Landslide Displacement Field

Landslides can have complex, spatially strongly inhomogeneous surface displacement fields with discontinuities from multiple active lobes that are deforming while failing on nested slip surfaces at different depths. For synthetic aperture radar interferometry (InSAR), particularly at lower resolutions, these characteristics can cause significant aliasing of the wrapped phase. In combination with steep terrain and seasonal snow cover, causing layover and temporal decorrelation, respectively, traditional phase unwrapping can become unfeasible, even after topographic phase contributions have been removed with an external high-resolution digital surface model (DSM). We present a novel method: warp demodulation that reduces the complexity of the phase unwrapping problem for noisy and/or aliased, low-resolution interferograms of discontinuous landslide displacement. The key input to our warp demodulation method is a single (or several) reference interferogram(s) from a high-resolution sensor mode such as TerraSAR-X Staring Spotlight with short temporal baseline and good coherence to allow localization of phase discontinuities and accurate unwrapping. The task of constructing suitable phase surfaces to approximate individual to-be-demodulated interferograms from the reference interferogram is made difficult by strong and spatially inhomogeneous temporal, seasonal, and interannual variations of the landslide with individual lobes accelerating or decelerating at different rates. This prevents using simple global scaling of the reference. Instead, our method uses an irregular grid of small patches straddling strong spatial gradients and phase discontinuities in the reference to find optimum local scaling factors that minimize the residual phase gradients across the discontinuities after demodulation. Next, for each to-be-demodulated interferogram, from these measurements we interpolate a spatially smooth global scaling function, which is then used to scale the (discontinuous) reference. Demodulation with the scaled reference leads to a residual phase that is also spatially smooth, allowing it to be unwrapped robustly after low-pass filtering. A key assumption of warp demodulation is that the locations of the phase discontinuities can be mapped in the reference and that they are stationary in time at the scale of the image resolution. We carry out extensive tests with simulated data to establish the accuracy, robustness, and limitations of the new method with respect to relevant parameters, such as decorrelation noise and aliasing along phase discontinuities. A realistic parameterization of the method is demonstrated for the example of the Fels Glacier Slide in Alaska using a recent late-summer high-resolution staring spotlight interferometric image pair from TerraSAR-X to demodulate. We show warp demodulation results for also recent but early-summer, partially incoherent interferograms of the same sensor, as well as for older and coarser aliased interferograms from RADARSAT-2, ALOS-1, and ERS.


Introduction
Interferometric synthetic aperture radar [1,2] (InSAR) processing techniques for measuring the Earth's surface deformation have evolved dramatically over the last decades.First generation techniques to extract the deformation phase, based on the simple removal of the topographic phase components in pairwise interferograms [3][4][5] have been complemented by second generation time series analysis techniques, which exploit temporal networks of interferograms for the simultaneous statistical modeling of the topographic, atmospheric, and deformation phase components, with [6][7][8][9] or without prior phase unwrapping [10][11][12].Lastly, corresponding third generation techniques have been developed for the wrapped [13] and unwrapped case [14,15], respectively, which use additionally spatially adaptive filtering to improve coherence, while preserving spatial resolution.Also, for the phase unwrapping of interferometric networks, more powerful spatio-temporal analogues [16] have been developed beyond the more basic pairwise algorithms [17].The problem of unwrapping the InSAR phase has also been simplified already via partial demodulation of the wrapped phase, both through the availability of higher resolution global digital elevation models (DEMs) (e.g., [18,19]) and the successful introduction of deformation models for specific types of phenomena causing surface displacement.Examples where additional deformation modeling and demodulation can lead to small phase residuals that ease phase unwrapping considerably include volcanoes, earthquakes, and glaciers [20,21].
In contrast, despite the technical sophistication in processing techniques and the availability of high-resolution DEMs, recovering from InSAR data the highly temporally variable and spatially discontinuous deformation fields of large complex landslides have remained challenging to both time series methods and phase unwrapping algorithms [22].The inability to achieve a small phase residual via modeling of the landslide deformation is mainly due to the lack of accurate geophysical models in combination with other compromising factors, such as poor or patchwise-variable phase quality from decorrelation (e.g., due to precipitation, seasonal snow melt, and vegetation) or from surface disintegration caused by the deformation itself, as well as a multitude of small layover and shadow areas caused by rugged, steep terrain.
We have developed a simple and robust phase demodulation method-warp demodulation-that uses high-resolution interferometric data to unwrap the phase of interferograms that are partially incoherent and/or aliased due to coarser spatial resolution.The method relies crucially on the availability of at least one high-resolution coherent template where active faults can be delineated.The parameterization of warp demodulation has been specifically tailored to the example of the Fels Glacier Slide, Alaska, and consequently, the applicability of the method is restricted currently to the diversity provided by this example.The deformation field of the Fels Glacier Slide has pronounced discontinuities caused by multiple active lobes failing on nested slip surfaces at various depths.There also is strong and spatially inhomogeneous temporal variation of the deformation field, as well as significant spatio-temporal coherence variations caused by precipitation and melting of seasonal snow cover.However, the vegetation cover of the Fels Glacier Slide is low (alpine tundra), and the terrain ruggedness is only moderate.To provide a wider context for the current version of the warp demodulation method, we summarize in Table 1 the general complicating factors to the successful phase unwrapping of interferograms of landslide motion and how far these are realized by our case study on the Fels Glacier Slide.We believe our method to be transferable with some adaptation to other landslide scenarios.
The scope of this paper is also limited to only the spatial demodulation of the discontinuous landslide displacement field of the Fels Glacier Slide.Extending our warp demodulation method to include a direct coupling with the strong temporal non-linearity originating from the pronounced temporal variations in the seasonal and inter-annual landslide drivers observed for our example and most other landslides is the subject of follow-up work.The presented demodulation solution can be regarded as a first step towards a robust and automated InSAR time series analysis solution optimized for complex landslide deformation.
The paper is organized as follows.Section 2 outlines the novel phase demodulation method.In Section 3, we first present results from a sensitivity analysis of the new method using simulated data.We then demonstrate the method for the Fels Glacier Slide using recent high-resolution Staring Spotlight [23,24] TerraSAR-X (TSX-ST) data to demodulate also recent but seasonally less coherent data from TerraSAR-X, as well as earlier, coarser resolution aliased data from RADARSAT-1, ALOS-1, and ERS.The results are discussed in terms of capabilities and limitations, as well as in the context of planned extensions of the new method in Section 4.

Method
Our initial idea to construct appropriate demodulation surfaces that reduce the complexity of the phase unwrapping problems for single interferograms of spatially discontinuous landslide motion was to use a wireframe model similar to [25].The wireframe model consists of an adjustable parametric surface constructed over a triangular mesh.The mesh boundary includes cuts at the active fault line locations of the slide, which enable the model to fit the displacement field first-order differentiable except at its discontinuities.In this study, we use a related but much simpler warp demodulation model that directly matches an unwrapped reference interferogram to the to-be-demodulated interferograms at discrete locations along the active faults.The locally optimized match factors are then interpolated to fit the reference to the wrapped interferograms with a global multiplicative warp function.Compared to the full wireframe model, our warp demodulation model has the key advantage of being straightforward to parametrize robustly.The disadvantage of the warp demodulation model of leaving a larger demodulation residual that can still be wrapped is largely neutralized by this residual being continuous and slowly varying spatially.As a result, the residual-in contrast to the original discontinuous interferogram-can be unwrapped easily and accurately after applying a simple boxcar low pass (multi-looking).The steps of our method are schematically depicted at the top of Figure 1 and are explained in more detail as follows: Step 1 Detect Active Fault Locations.Active fault locations can be mapped robustly and accurately by using a curvature criterion (second-order derivative; minimum threshold on the largest eigenvalue and contrast of eigenvalues of the two-dimensional (2D) Hessian matrix of factor >10) applied to a smoothed version of the reference interferogram.This, or a similar discontinuity detection approach, must be used by the wireframe model, which requires continuous parameterization of each active fault line in its entirety.For the warp demodulation model, it is sufficient for each fault line to be represented by a (generally coarser and irregularly spaced) representative set of points.Such a set can be found by subsampling the output from the described curvature criterion discontinuity detector.This approach was used for the Fels Glacier Slide example (see Section 3.2).Alternatively, it is possible to digitize the point set representing the faults manually on the reference interferogram (as was done for the sensitivity analysis in Section 3.1); this approach may be generally more robust, particularly if a less coherent and high-resolution reference is used.
Step 2 Find Local Match Factors.The reference is matched to the to-be-demodulated interferograms within circular patches straddling the active faults at the pointset locations found in the first step.The radius of the patches is an important parameter of the model whose influence on the accuracy and robustness of the demodulation is discussed later.Optimum multiplicative factors to scale the reference before demodulation are obtained now for each patch by maximizing the (signal-biased) coherence estimate of the residual inside the patch.While this type of match (coherence optimization) essentially minimizes the dynamic range of the coherent phase residual of the patch (incoherent phase noise inside the patch does not change during the optimization and thus does not affect it) it is insensitive to the phase offset between the interferogram and the scaled reference.This is a key feature of our warp demodulation method, which favors spatial smoothness over smallness of the demodulation residual.
Step 3 Thin Plate Spline Fit.The vector of optimum match factors at the point set locations is now interpolated to form a continuous function of the interferogram coordinates using a thin plate spline fit.This function is then multiplied with the reference to warp the latter into an optimum demodulation phase surface for the interferogram.This demodulation phase is then wrapped, and its complex conjugate is multiplied with the to-be-demodulated interferogram to form the residual phase.
Step 4 Unwrap Residual.The residual phase, optimally smooth and slowly varying spatially, is multi-looked and unwrapped using standard Minimum Cost Flow (MCF) techniques [17].
The unwrapped residual is added to the constructed (unwrapped) demodulation phase of step 3 to form the final unwrapped interferogram.

Results
In this section, we present results from a sensitivity analysis of the new method, as well as demonstration results using actual data.Table 2 contains a list of symbols and acronyms used in Section 3. We want to point out explicitly that our warp demodulation method assumes the existence of a reference interferogram of "good quality", quantified as high coherence and unobscured by layover with an MCF unwrapping result that is correct except for perhaps localized small errors (restricted to the vicinity of a few pixels around minor incoherent areas).Therefore, within the scope of the present study, the effects arising from shortcomings of the reference have not been included in the sensitivity analysis.The warp demodulation method has several unique characteristics: (i) as already mentioned, the method is designed to produce a smooth demodulation residual; the mean phase offset for each patch between the reference and the to-be-demodulated interferogram deliberately is not considered in finding the optimum local matching factors; (ii) phase structures inside the patch instead are used for the optimization; primarily, the phase discontinuity along the fault line dissecting the patch is used, but secondarily, the gradients and curvature features of the displacement phase of the landslide lobes are also used.The ability to use this secondary information allows the method to succeed in finding the correct demodulation surface, even if phase aliasing occurs across the fault due to the acceleration of a lobe along its downslope axis.The chosen patch radius plays an important role in how robustly the method delivers on this ability, and this is discussed in Section 3.1, where we present a full parametric sensitivity analysis of the method; (iii) our method relies on approximate stationarity of the phase discontinuities between the reference and the interferogram: significant shifting of active faults leads to artifacts also discussed in Section 3.1.The method biases the local scale factor to zero for a "flat" coherence optimum, which handles the case correctly that when an active fault in the reference is inactive in the interferogram, the resulting residual is smooth.However, for the opposite case where a fault in the interferogram is not present in the reference, the coherence optimization will be insensitive and cannot diminish the signal (residual will remain discontinuous).The wireframe model has the same problem; if new active faults are generated, then the reference needs to be updated to include the new faults as cuts.As there, (iv) the warp demodulation method can be generalized straightforwardly by including several references appropriately spaced temporally and selecting the one close enough in time for each to-be-demodulated interferogram to capture adequately the formation (and cessation) of active faults.

Results
In this section, we present results from a sensitivity analysis of the new method, as well as demonstration results using actual data.Table 2 contains a list of symbols and acronyms used in Section 3. We want to point out explicitly that our warp demodulation method assumes the existence of a reference interferogram of "good quality", quantified as high coherence and unobscured by layover with an MCF unwrapping result that is correct except for perhaps localized small errors (restricted to the vicinity of a few pixels around minor incoherent areas).Therefore, within the scope of the present study, the effects arising from shortcomings of the reference have not been included in the sensitivity analysis.

Sensitivity Analysis with Simulated Data
The warp demodulation method is applied to a dataset of simulated InSAR observations for two different landslide scenarios:

•
Single-lobe: the deforming active lobe is modeled by a semi-ellipse (vertical semi-axis of 80 pixels, horizontal axis of 100 pixels, Figure 2a).

•
Multi-lobe: three nested semi-elliptical lobes are simulated (see Figure 2b for dimensions).The goal of this analysis is to assess the robustness and limitations of warp demodulation with respect to the following landslide deformation and/or demodulation parameters (see also Figure 3):     to zero if ∆ < 1 fringe/lobe.Above this threshold, the RMSE tends to increase with ∆ until the latter reaches approximately 1.6-1.8fringes/lobe, then the two quantities become again uncorrelated.Such an increase is markedly abrupt for f 1 fringe/lobe (from 0 to 0.96 rad as ∆ increases from 1 to 1.5 fringes/lobe).
For ∆ values bigger than 1 fringe/lobe, the RMSE is also negatively associated with f .For instance, RMSE at ∆ 1.5 fringes/lobe is seen to increase by 0.85 rad in response to a f reduction by 50%.
Logistic gradient: If f 1 fringes/lobe, the RMSE remains approximately equal to zero until approximately ∆ = 1 fringe/lobe, when it abruptly increases, reaching a maximum of 1.05 rad at ∆ = 1.7 fringes/lobe.The same behavior with ∆ , although with a less abrupt increase, is also found for f 1.5 fringes/lobe.If 2 f 4 fringes/lobe, the RMSE shows a positive association with ∆ (Pearson's correlation coefficient 0.93 ); in this case, the statistics appear substantially insensitive to f and range from a minimum of 0.2 rad (at ∆ = 0.3 fringes/lobe) to a maximum of 0.5 rad (at ∆ = 1.9 fringes/lobe).

Effect of Atmosphere-Like Disturbances
We simulated a two-dimensional atmospheric phase error surface (GAU) that shares the spatial characteristics of the dynamic atmosphere.GAU is assumed to be zero-mean Gaussian-distributed (with standard deviation σ) and is smoothed via an 11-pixel uniform filter.The to-be-demodulated interferogram IFG' now has the form: In Tables 4 and 5, we reported the mean RMSE of the demodulated unwrapped phase computed for a population of 50 randomly generated IFG' observations at three different values.The statistics are obtained by assuming either a ramp or a logistic ( 4) along-lobe deformation gradient.
The mean RMSE shows a trend with ∆ that is in agreement with the results of Section 3.1.1(constant until ∆ = 1 fringe/lobe, then abruptly increasing).Compared to Section 3.1.1,all RMSE Ramp gradient: The RMSE remains substantially insensitive to f REF and is approximately equal to zero if ∆ f < 1 fringe/lobe.Above this threshold, the RMSE tends to increase with ∆ f until the latter reaches approximately 1.6-1.8fringes/lobe, then the two quantities become again uncorrelated.Such an increase is markedly abrupt for f REF = 1 fringe/lobe (from 0 to 0.96 rad as ∆ f increases from 1 to 1.5 fringes/lobe).
For ∆ f values bigger than 1 fringe/lobe, the RMSE is also negatively associated with f REF .For instance, RMSE at ∆ f = 1.5 fringes/lobe is seen to increase by 0.85 rad in response to a f REF reduction by 50%.
Logistic gradient: If f REF = 1 fringes/lobe, the RMSE remains approximately equal to zero until approximately ∆ f = 1 fringe/lobe, when it abruptly increases, reaching a maximum of 1.05 rad at ∆ f = 1.7 fringes/lobe.The same behavior with ∆ f , although with a less abrupt increase, is also found for f REF = 1.5 fringes/lobe.If 2 < f REF < 4 fringes/lobe, the RMSE shows a positive association with ∆ f (Pearson's correlation coefficient R = 0.93); in this case, the statistics appear substantially insensitive to f REF and range from a minimum of 0.2 rad (at ∆ f = 0.3 fringes/lobe) to a maximum of 0.5 rad (at ∆ f = 1.9 fringes/lobe).

Effect of Atmosphere-Like Disturbances
We simulated a two-dimensional atmospheric phase error surface (GAU) that shares the spatial characteristics of the dynamic atmosphere.GAU is assumed to be zero-mean Gaussian-distributed (with standard deviation σ) and is smoothed via an 11-pixel uniform filter.The to-be-demodulated interferogram IFG' now has the form: In Tables 4 and 5, we reported the mean RMSE of the demodulated unwrapped phase computed for a population of 50 randomly generated IFG' observations at three different σ values.The statistics are obtained by assuming either a ramp or a logistic (κ = 4) along-lobe deformation gradient.The mean RMSE shows a trend with ∆ f that is in agreement with the results of Section 3.1.1(constant until ∆ f = 1 fringe/lobe, then abruptly increasing).Compared to Section 3.1.1,all RMSE measurements of Tables 4 and 5 are now affected by a positive offset that is positively associated with σ (e.g., compared to the case where σ = 3, the RMSE at σ = 11 is on average 0.44 rad greater).

Impact of the Demodulation Patch Radius
To assess the performance of the warp demodulation approach with respect to the radius of the demodulation patches r, the RMSE of the demodulated unwrapped phase is plotted against r when either a ramp or a logistic deformation gradient is imposed (Figure 5a,b, respectively).All simulations are performed using six patches located on the main lobe discontinuity.Logistic gradient (∆ f = 0.75, f REF = 4): The RMSE is consistently smaller at κ = 2 than κ = 4 and is seen to decrease with increasing r when the latter is below 0.2 (κ = 2) or 0.35 (κ = 4).Above these thresholds, the RMSE remains consistently below 0.1 rad and positively associated with r (R = 0.84 at κ = 2; R = 0.96 at κ = 4).
The warp demodulation is also applied to a multi-lobe scenario with ∆ f = 0.75 fringes/lobe and f REF = 1.5 fringes/lobe (ramp gradient).As shown in Figure 6, a small radius r = 0.1 results in localized phase errors within the main lobe (RMSE = 0.206 rad).When 0.3 ≤ r ≤ 0.5, the phase gradient along the main lobe is accurately reconstructed.However, a number of phase artifacts within the secondary lobes can now be observed.The number of artifacts, and hence the RMSE, appear to increase with increasing r (r = 0.3 : RMSE = 0.181 rad, r = 0.5 : RMSE = 1.318 rad).

Real Data Example: Fels Glacier Slide
In the following, we give a demonstration of the warp demodulation method for a representative large, complex landslide test site, the Fels Glacier Slide in Alaska.The Fels Glacier Slide [27] occupies an instable, south-facing slope above and immediately adjacent to the lower Fels Glacier.The Fels Slide is about 3 km from the Alaska Pipeline corridor to the west and about 2 km from the Denali Fault to the south; earthquakes generated by the Denali Fault, such as the M7.9 one in November 2002, significantly affect the slide dynamics.Other obvious drivers of the dynamics of the slide are meteorology (mostly snow melt and rain) and the ongoing recession of the Fels Glacier, which additionally destabilizes the toe of the slide.On the surface, the deforming area of the slide is about 3 km wide and has approximate fall line and elevation extents of 2 km and 600 m, respectively.The mechanism of the deformation is a deep-seated (depth and shape of the failure plane unknown) gravitational mass deformation forming two large lobes separated east-west by a central dissecting gully.Overlain onto these two deep-seated main lobes are several smaller more surficial deforming lobes that cause strong discontinuities in the deformation field where the surface ruptures along extensive downward curving arcuate fissures.While the deformation is classified as creep-type, movement amounts can easily exceed 30 cm per year for the fastest areas of the slide.
As inputs to our demodulation method, to construct the high-resolution reference "model" we use a TSX-ST image pair 3-14 August 2015 (0.6 m spatial resolution; shown in Figure 1), and for the

Real Data Example: Fels Glacier Slide
In the following, we give a demonstration of the warp demodulation method for a representative large, complex landslide test site, the Fels Glacier Slide in Alaska.The Fels Glacier Slide [27] occupies an instable, south-facing slope above and immediately adjacent to the lower Fels Glacier.The Fels Slide is about 3 km from the Alaska Pipeline corridor to the west and about 2 km from the Denali Fault to the south; earthquakes generated by the Denali Fault, such as the M7.9 one in November 2002, significantly affect the slide dynamics.Other obvious drivers of the dynamics of the slide are meteorology (mostly snow melt and rain) and the ongoing recession of the Fels Glacier, which additionally destabilizes the toe of the slide.On the surface, the deforming area of the slide is about 3 km wide and has approximate fall line and elevation extents of 2 km and 600 m, respectively.The mechanism of the deformation is a deep-seated (depth and shape of the failure plane unknown) gravitational mass deformation forming two large lobes separated east-west by a central dissecting gully.Overlain onto these two deep-seated main lobes are several smaller more surficial deforming lobes that cause strong discontinuities in the deformation field where the surface ruptures along extensive downward curving arcuate fissures.While the deformation is classified as creep-type, movement amounts can easily exceed 30 cm per year for the fastest areas of the slide.
As inputs to our demodulation method, to construct the high-resolution reference "model" we use a TSX-ST image pair 3-14 August 2015 (0.6 m spatial resolution; shown in Figure 1), and for the to-be-demodulated lower resolution interferograms, we have selected three representative examples, one each from ERS, ALOS-1, and RADARSAT-2 capturing different resolution (30 m, 15 m, and 3 m, respectively), seasonal and annual separation, coherence conditions, and interferogram periods with respect to the reference (Table 6).The TSX-ST reference is unwrapped using MCF and re-projected to the line-of-sights of the respective other sensors.The re-projection assumes an approximate motion direction following the average fall line as derived from a spatially smoothed version (500 m scale) of the external DSM used for the topographic corrections and is generally robust due to the similarity of acquisition geometries of the different sensors.Track angles in Table 6 are very similar, as expected for space-borne data of the same orbit direction; look angles are more different but still only vary by less than 20 degrees, which for only moderately rugged topography (as is the case for the Fels Glacier Slide) does not lead to stark layover/shadow contrasts for the different incidence angles, guaranteeing a smooth re-projection map.After re-projection of the template, the warp demodulation method is then applied as introduced in Section 2 and analyzed in Section 3.1.Standard MCF unwrapping (using coherence costs and masking areas with coherence below 0.3) results of four selected test interferograms over the Fels Glacier Slide are shown in Figure 7, subfigures (a)-(d).For easy comparison, subfigures (e)-(h) of Figure 7 show the final corresponding warp demodulation results presented in detail in Figures 8-11 (lower right panels).Problems from erroneous branch cuts made by the MCF algorithm are immediately apparent for both the ERS case (within the main lobe, Figure 7a) and most strikingly for the TSX case (Figure 7d).For the ALOS-1 and RSAT-2 cases (Figure 7b,c, respectively), MCF performs better; however, more localized branchcut problems are still easily identifiable.
Remote Sens. 2018, 10, x FOR PEER REVIEW 12 of 17 Also, for the other two cases (Figures 10 and 11) involving sensors with considerably coarser resolution, our method still performs satisfactorily in producing phase residuals that are both spatially smooth and sufficiently small to reconstruct the unwrapped phase of the test interferograms.Reconstruction is not perfect, firstly due to the error sources identified in Section 3.1 (e.g., migration and introduction of deformation discontinuities), plus a few additional ones (such as small unwrapping errors in the reference in some small layover regions but particularly in the fastmoving area where the eastern lobe meets the Fels Glacier).Errors in the reconstruction are by and large limited to this fastest moving part of the slide (most visible in Figure 10e, the unwrapped residual of the ALOS-1 case), where the reference likely would also not be representative due to separation in time even in the absence of any phase unwrapping errors.

Discussion
The sensitivity analysis of the warp demodulation method in Section 3.1 using synthetic data identifies a number of error sources assuming a perfect reference.Our method, for obvious reasons, cannot reconstruct areas where the reference does not exist.However, in areas where the highresolution interferogram used to generate the demodulation reference is of high quality and has been correctly unwrapped, our method is robust with respect to almost all potential error sources; the ones analyzed in Section 3.1 and others associated with the described processing steps, such as the reprojection of the unwrapped reference to the line-of-sight of the to-be-demodulated interferogram.Of particular interest is the ability of our method to recover the unwrapped phase correctly across a The choice of demodulation patch sizes reflects our reasoning on the different resolutions and atmospheric and random noise phase error levels of the different sensors.For example, ERS has the largest radius (350 pixels); ERS has comparable atmospheric error to RADARSAT-2, but its random noise is higher due to factor 10 lower resolution.ALOS-1 (patch radius 250 pixels), covering one entire year, has random noise only slightly lower than ERS over one month, but tropospheric atmospheric noise levels in L-band are expected to be much lower than for ERS.Due to its considerably high-resolution, the RADARSAT-2 patch radius is the smallest (200 pixels).
The findings from the real data experiments indicate that, unlike standard MCF unwrapping (Figure 7d), the warp demodulation method performs creditably over the TSX-ST interferogram of May 2016 (Figure 8) despite the extensive (around 50 percent of area) low coherence areas of the original interferogram.Further, the findings confirm (most obvious in the 48-day C-band interferogram of Figure 9) that our method can handle significant aliasing in IFG that "hides" active discontinuities through a near integer 2-pi phase ambiguity; warp demodulation still finds the correct FAC values across these discontinuities by matching between IFG and REF the spatial texture of the phase surfaces on both sides of the fault.RES in Figure 9d is smooth and within a single fringe, which indicates a successful warp demodulation.Also, for the other two cases (Figures 10 and 11) involving sensors with considerably coarser resolution, our method still performs satisfactorily in producing phase residuals that are both spatially smooth and sufficiently small to reconstruct the unwrapped phase of the test interferograms.Reconstruction is not perfect, firstly due to the error sources identified in Section 3.1 (e.g., migration and introduction of deformation discontinuities), plus a few additional ones (such as small unwrapping errors in the reference in some small layover regions but particularly in the fast-moving area where the eastern lobe meets the Fels Glacier).Errors in the reconstruction are by and large limited to this fastest moving part of the slide (most visible in Figure 10e, the unwrapped residual of the ALOS-1 case), where the reference likely would also not be representative due to separation in time even in the absence of any phase unwrapping errors.

Discussion
The sensitivity analysis of the warp demodulation method in Section 3.1 using synthetic data identifies a number of error sources assuming a perfect reference.Our method, for obvious reasons, cannot reconstruct areas where the reference does not exist.However, in areas where the high-resolution interferogram used to generate the demodulation reference is of high quality and has been correctly unwrapped, our method is robust with respect to almost all potential error sources; the ones analyzed in Section 3.1 and others associated with the described processing steps, such as the re-projection of the unwrapped reference to the line-of-sight of the to-be-demodulated interferogram.Of particular interest is the ability of our method to recover the unwrapped phase correctly across a lobe discontinuity, even if the phase has several wraps when followed along the discontinuity delineating the lobe.This ability depends on the presence of a significantly strong spatial texture (local gradient or curvature) in the phase that can be matched between reference and to-be-demodulated interferogram to break the two-pi ambiguity of the phase jump across the discontinuity.To match spatial textures robustly, a large enough patch size is critical.However, a very large size can also cause patches to overlap with multiple discontinuities, where they are closely spaced, which prevents a smooth residual should differential deformation changes occur for the different discontinuities with respect to the reference.The analysis in Section 3.1.3revealed quantitatively how patch size trades off quantitatively for realistic lobe configurations.Our warp demodulation method is generally robust for a range of patch sizes between a half and a fifth of the main lobe diameter.An enhanced future version of our method will use an adaptive patch size that takes the variable spacing of phase discontinuities (slide faults) into account to combine the increase in robustness with patch size with preserved high-resolution where several discontinuities are nested close together.
A noteworthy exception to the generally robust behavior of the warp demodulation method is when the reference gets "out of date" (i.e., the number of discontinuities of the deformation field identified in the reference does not match that in the to-be-demodulated interferogram (this can be considered a more extreme case of the subtle shifting of discontinuities included in the analysis of Section 3.1)).For the case of a qualitative mismatch of discontinuities, where either (i) an extra discontinuity or (ii) a discontinuity identified in the reference is not active in the interferogram, there is a pronounced asymmetry.For case (ii), the warp factor will down-scale the reference locally at the inactive discontinuity; the warp modulation method is robust to this case (=smooth residual).However, for case (i), the residual will be insensitive to the warp factor and the residual will be locally discontinuous.Thus, if too many discontinuities are missing in the reference and there are additional high-resolution datasets available closer in time, then a new reference should be constructed and used instead.
As mentioned before, our method is easily generalized to several references and can also be used iteratively (i.e., successfully warp-demodulated interferograms can themselves be used as reference).A full generalization of our method to a wireframe representation in space and time is considered as a next step.For a slide with a deep, historic InSAR dataset, this would then include an explicit parameterization of the spatio-temporal pattern of the active faults of the landslide, including formation, migration, and healing of faults.The patch measurements carried out with our original warp demodulation concept presented in Section 2 could be used to make the explicit fault parametrization converge.Once found, the explicit fault parametrization can then be used to demodulate/unwrap the full time series of interferograms (as discontinuous aliased phase surfaces) in a single optimization step (as per the original wireframe concept).
Besides the robustness of the method, the interpretation of the two results per interferogram (i) warp factor map and (ii) residual phase map is also of interest for discussion.As our method targets smoothness but not necessarily the overall smallest size of the phase residual after demodulation, neither the residual, nor the warp factor maps can be easily interpreted quantitatively as a measure of bulk slide activity relative to that of the reference.Qualitatively, a smooth residual shows persistence of the fault pattern between the interferogram and the reference.Also, a factor map in the 1-range plus a small residual indicates that the deformation field of interferogram and the reference are very similar.However, beyond that-despite obviously containing important information about the (temporal) changes of the slide dynamics-the interpretation of the pair residual/warp factor map is complex.The studying of a variety of real data examples, such as the ones in Figures 8-11, has revealed that the warp factor difference forms one measured change in the shallow (discontinues) deformation field of the slide, while the size of the residual indicates a change in the deep deformation pattern between the interferogram and the reference.More experiments with simulated data are planned to understand this connection more quantitatively.

Conclusions
We have developed a new demodulation method (warp demodulation) capable of producing smooth phase residuals to improve the phase unwrapping performance for older, lower resolution, aliased interferograms of complex, multi-lobed, discontinuous landslide deformation.Warp demodulation hinges on the identification of the discontinuity (active fault) pattern on the slide

Figure 1 .
Figure 1.Top: the four steps of the warp demodulation method; IFG: to-be-demodulated wrapped interferogram; REF: high-resolution unwrapped reference interferogram; FAC: interpolated factor map; and RES: wrapped residual.Bottom: real data example of reference interferogram (TerraSAR-X Staring Spotlight, 0.6 m ground resolution), Fels Glacier Slide, Alaska.

Figure 1 .
Figure 1.Top: the four steps of the warp demodulation method; IFG: to-be-demodulated wrapped interferogram; REF: high-resolution unwrapped reference interferogram; FAC: interpolated factor map; and RES: wrapped residual.Bottom: real data example of reference interferogram (TerraSAR-X Staring Spotlight, 0.6 m ground resolution), Fels Glacier Slide, Alaska.

Figure 2 .
Figure 2. Sensitivity analysis with simulated data.Schematic depiction of the two simulated landslide scenarios: (a) single-lobe (semi-ellipse) and (b) multi-lobe (three nested semi-ellipses).The red dots identify the location of the demodulation patches (of radius ).

Figure 3 .
Figure 3. Sensitivity analysis with simulated data.Schematic depiction of the simulation experiments.SHIFT: robustness of the warp demodulation method against downward shifts Δh of the IFG active lobe; GRADIENT: algorithm sensitivity to along-lobe deformation differences ∆ between REF and IFG for ramp/logistic deformation gradients; GAUSSIAN: IFG affected by Gaussian noise; SPECKLE: IFG affected by fully developed speckle; RADIUS: Impact of patch radius on demodulation performance.

Figure 2 .
Figure 2. Sensitivity analysis with simulated data.Schematic depiction of the two simulated landslide scenarios: (a) single-lobe (semi-ellipse) and (b) multi-lobe (three nested semi-ellipses).The red dots identify the location of the demodulation patches (of radius r).

•Figure 2 .
Figure 2. Sensitivity analysis with simulated data.Schematic depiction of the two simulated landslide scenarios: (a) single-lobe (semi-ellipse) and (b) multi-lobe (three nested semi-ellipses).The red dots identify the location of the demodulation patches (of radius ).

Figure 3 .
Figure 3. Sensitivity analysis with simulated data.Schematic depiction of the simulation experiments.SHIFT: robustness of the warp demodulation method against downward shifts Δh of the IFG active lobe; GRADIENT: algorithm sensitivity to along-lobe deformation differences ∆ between REF and IFG for ramp/logistic deformation gradients; GAUSSIAN: IFG affected by Gaussian noise; SPECKLE: IFG affected by fully developed speckle; RADIUS: Impact of patch radius on demodulation performance.

Figure 4 .
Figure 4. Sensitivity analysis with simulated data, algorithm robustness against along-lobe deformation differences between REF and IFG.RMSE of the demodulated unwrapped phase as a function of the difference in the along-lobe deformation rate ∆ f f for four different REFs along-lobe deformation rates f .The statistics are obtained by assuming (a) a ramp along-lobe deformation gradient or (b) a logistic gradient (curve steepness 4).

Figure 4 .
Figure 4. Sensitivity analysis with simulated data, algorithm robustness against along-lobe deformation differences between REF and IFG.RMSE of the demodulated unwrapped phase as a function of the difference in the along-lobe deformation rate ∆ f = f IFG − f REF for four different REFs along-lobe deformation rates f REF .The statistics are obtained by assuming (a) a ramp along-lobe deformation gradient or (b) a logistic gradient (curve steepness κ = 4).

Figure 6 .
Figure 6.Sensitivity analysis with simulated data, impact of the demodulation patch radius on the algorithm performance (multi-lobe scenario).IFG: To-be-demodulated interferogram (wrapped); REF: reference interferogram (unwrapped); FAC: demodulation factor; RES: demodulation residual (unwrapped); UNW: demodulated interferogram.FAC, RES, and UNW are obtained for three different values.The phase errors in UNW are highlighted with white dashed lines.

Figure 6 .
Figure 6.Sensitivity analysis with simulated data, impact of the demodulation patch radius r on the algorithm performance (multi-lobe scenario).IFG: To-be-demodulated interferogram (wrapped); REF: reference interferogram (unwrapped); FAC: demodulation factor; RES: demodulation residual (unwrapped); UNW: demodulated interferogram.FAC, RES, and UNW are obtained for three different r values.The phase errors in UNW are highlighted with white dashed lines.

Figure 11 .
Figure 11.Warp demodulation result for a representative ERS interferogram (C-Band, 30 m ground range resolution, 35-day repeat cycle) covering 1 cycle.(a) IFG-ERS, to-be-demodulated interferogram (wrapped); (b) REF-TSX, reference interferogram (unwrapped); (c) FAC, demodulation factor; (d) RES: demodulation residual (wrapped); (e) RES-UNW: demodulation residual (unwrapped); and (f) UNW-ERS, demodulated interferogram.The unwrapped phase in (b,e,f) is clipped to ±2π.The demodulation patch radius is 350 pixels.Respective warp demodulation results for the same test sites are presented in Figures 8-11 (as just stated the unwrapped interferogram (UNW) of the lower right sub-figures are also shown in Figure 7e-h).The test case of Figure 8 (TSX-ST 2016) selects a partially incoherent interferogram with strong seasonal deformation contrast but with the same resolution and close in absolute time to the

Table 1 .
Complicating factors to the successful phase unwrapping/demodulation of landslide interferograms and the relevance of the Fels Glacier Slide example.

Table 2 .
List of parameters used in Section 3.
REF Along-lobe deformation rate in the REF interferogram; units: fringes/lobe f INF Along-lobe deformation rate in the INF interferogram.units: fringes/lobe ∆ f Difference (between REF and IFG) in the along-lobe deformation rate; units: fringes/lobe r Demodulation patch radius; units: fraction of the main lobe's diameter

Table 3 .
Sensitivity analysis with simulated data.Overview of the simulation experiments and results.

Table 4 .
Sensitivity analysis with simulated data, effect of atmosphere-like disturbances.RMSE (in radians) of the demodulated unwrapped phase at three different σ values (ramp gradient).

Table 5 .
Sensitivity analysis with simulated data, effect of atmosphere-like disturbances.RMSE (in radians) of the demodulated unwrapped phase at three different σ values (logistic gradient).

Table 6 .
Synthetic aperture radar (SAR) datasets used to test the warp demodulation method.

Table 6 .
Synthetic aperture radar (SAR) datasets used to test the warp demodulation method.