Satellite interferometry for mapping surface deformation time series in one, two and three dimensions: A new method illustrated on a slow-moving landslide

Satellite interferometry for mapping surface deformation time series in one, two and three dimensions: A new method illustrated on a slow-moving landslide

Introduction

Surface deformation information is a key parameter for characterization of natural and anthropogenic processes such as earthquakes, volcanic eruptions, landslides, glacier flows, and subsidence and uplift due to resource exploitation (e.g. Colesanti and Wasowski, 2006; Wasowski and Bovenga, 2014; Elliott et al., 2016; Bayer et al., 2017; Dong et al., 2018; Li et al., 2019). To this end, space-borne Differential Interferometric Synthetic Aperture Radar (DInSAR) allows the derivation of surface deformation with wide spatial coverage, millimetre precision, high temporal and spatial resolutions and all-day and all-weather working capabilities (Massonnet and Feigl, 1998; Rosen et al., 2000).

However, conventional DInSAR technique measures only one component of the surface deformation – in the satellite’s line-of-sight (LOS) (Wright et al., 2004; Hu et al., 2014). This makes interpretation of DInSAR measurements challenging and narrows the understanding of the mechanisms and dynamics of the deformation processes at work. Different methods were therefore developed to provide direct measurements of the two (2D, east and vertical) (e.g. Samsonov and d’Oreye, 2012, Samsonov and d’Oreye, 2017) or the three (3D, north, east and vertical) components of surface deformation at one point in time (e.g. Wright et al., 2004; Gourmelen et al., 2011; Ng et al., 2012; Raucoules et al., 2013; Shi et al., 2018). The most straightforward strategy for deriving 3D surface deformation is to combine multiple DInSAR LOS measurements with different viewing geometries. However, trying to derive 3D deformation from multi-orbit satellite SAR observations alone poorly constrains the north component of motion (Fuhrmann and Garthwaite, 2019). Indeed, due to the small inter-satellite variation in SAR satellite orbits (i.e. near-polar orbit), this approach may be applied solely for a very specific case – in high latitude areas where the azimuth angle diversity is higher (e.g. Gray, 2011; Hu et al., 2014). Larger diversity in acquisition geometry can be obtained using airborne DInSAR, that allows acquisition with any desired azimuthal direction (e.g. Delbridge et al., 2016; Handwerger et al., 2019). However, the cost of such acquisitions limits the temporal sampling. As a results of these difficulties, the retrieval of 3D surface deformation therefore commonly consists of merging conventional DInSAR LOS measurements with other types of measurements, that may rely either on SAR imagery or ancillary data, such as GNSS measurements (e.g. Samsonov and Tiampo, 2006; Samsonov et al., 2007, Samsonov et al., 2008; Hu et al., 2014).

East and vertical deformation components are often derived from conventional DInSAR applied to two (or more) different satellite tracks, that are then combined with outputs from either speckle tracking Wright et al., 2004; De Michele et al., 2010; Raucoules et al., 2013; Shi et al., 2018) or Multiple Aperture Interferometry (MAI) (e.g. Gourmelen et al., 2011) to retrieve the north component of deformation (Hu et al., 2014). With the speckle tracking technique, range and azimuth deformation are calculated from the offset of a given point on the ground using local correlation analysis between two images (Michel et al., 1999; Leprince et al., 2007). This method is applied to either the amplitude (Michel et al., 1999; Pathier et al., 2006) or the phase information (Gray et al., 2001; Pattyn and Derauw, 2002) of the SAR images. While it is not affected by signal saturation with high deformation gradient, speckle tracking precision is much lower, with values between 1/30th to 1/10th of the image pixel size, than the millimetre-scale precision that can be achieved with DInSAR (Leprince et al., 2007; Raucoules et al., 2013; Hu et al., 2014). This implies that the use of the speckle tracking technique is often limited to the case of relatively large deformation. The other method, MAI, aims at extending the applicability of DInSAR to along-track deformation based on split-beam DInSAR processing (Bechor and Zebker, 2006; Jung et al., 2009; Gourmelen et al., 2011). Bechor and Zebker (2006) found that for interferograms with coherence of 0.4 and higher, MAI outperforms speckle tracking analysis in precision. However, such high coherence is relatively uncommon, particularly in vegetated regions. Pepe et al. (2016); Pepe and Calo (2017) proposed to derive 3D deformation from ascending and descending data by applying higher order regularizations; this technique was used for mapping 3D deformation time series at the Big Island of Hawaii.

Another strategy used to overcome the constraint originating from the geometry of SAR observation is to make assumptions about deformation processes being studied (Hu et al., 2014). In the case of landslides and glacier flows, a reasonable assumption is that motion occurs parallel to the surface, although not necessarily downslope. By reducing the degrees of freedom from three to two, with this assumption one may resolve the three components of motion from (at least) one ascending and one descending DInSAR datasets. This approach has been demonstrated on glaciers in polar regions using single or few ascending and descending interferograms from the ERS-1/2 tandem mission (e.g. Joughin et al., 1998; Mohr et al., 1998; Gourmelen et al., 2011; Kumar et al., 2011). Notti et al. (2014); Kalia (2018); Hu et al., 2018, Hu et al., 2019 applied a similar principle on landslides, here projecting measured LOS velocity along the steepest downslope (one step further from assuming a movement parallel to the surface) to obtain estimates of the landslide motion. Assumption of motion along the steepest downslope may however be seen as too strict, as it does not allow consideration of lateral motion.

Time series are essential for studying surface deformation processes, which ideally calls for an integrated multi-sensor, multi-track and multi-temporal DInSAR measurements approach (Hu et al., 2013, Hu et al., 2014). However, all the previous methods discussed above, derived 3D surface deformation from either individual interferograms, or separate or interpolated ascending and descending datasets. Because of the mismatch in the acquisition times of ascending and descending datasets, these methods cannot be used to merge all available DInSAR data for producing a single set of temporally dense, 3D surface deformation time series (Samsonov and d’Oreye, 2012).

Here, we propose a novel version of the Multidimensional Small Baseline Subset (MSBAS) method (Samsonov and d’Oreye, 2012, Samsonov and d’Oreye, 2017), MSBAS-3D, that is able to produce 3D (north, east, and vertical) surface deformation time series from multi-track and multi-sensor DInSAR data. This new method aims at measuring surface deformation for processes involving motion parallel to the surface, such as landslides and glacier flows. The ability of MSBAS-3D to capture the full 3D deformation pattern and the benefits for the understanding of the processes at work are illustrated for a large, slow-moving, deep-seated landslide for which long DInSAR and dGNSS time series, as well as ground truth data, are available.

Methodology

The Multidimensional Small Baseline Subset (MSBAS) methodology is developed for simultaneous post processing of the hundreds to thousands of individual interferograms produced by a conventional DInSAR processing (Samsonov and d’Oreye, 2012, Samsonov and d’Oreye, 2017). Depending on the number of datasets, their acquisition parameters and specified control flags, the technique produces line-of-sight (MSBAS-1D), east and vertical (MSBAS-2D), and now north, east, and vertical (MSBAS-3D) surface deformation time series. The output of processing consists of cumulative time series for each acquisition epoch and linear deformation rate for each pixel that remains coherent throughout the entire time span. A brief derivation of all three techniques is provided below, with an emphasis on the new MSBAS-3D technique.

MSBAS-1D

In case of a single DInSAR dataset computed from data acquired by a sensor with an azimuth angle θ and an incidence angle ϕ, the LOS surface deformation time series is computed by applying the Small Baseline Subset (SBAS) technique (Berardino et al., 2002; Usai, 2003). The SBAS technique implemented in MSBAS-1D software is called when one or more DInSAR datasets with similar incidence and azimuth angles are provided. In matrix form MSBAS-1D can be written as

where A is a matrix constructed from the time intervals between consecutive SAR acquisitions, Vlos is a vector of the unknown LOS deformation velocities, Φ is a vector of observed DInSAR values, A+ is a pseudo-inverse of matrix A found by applying the Singular Value Decomposition (SVD), and dlosi is a cumulative LOS deformation at the time ti.

The LOS deformation velocity Vlos is related to the three components of the deformation velocity vector V with components VN, VE, VV through the following equation

where S = {SN, SE, SV} = {sinθ sin ϕ, − cos θ sin ϕ, cosϕ} is a LOS unit vector with north, east and vertical components SN, SE, SV.

MSBAS-2D

When both, ascending and descending DInSAR datasets acquired over the same area and time period are available, the MSBAS-2D technique computes east and vertical components of surface deformation time series. This solution is approximate because the contribution of the third, north component is neglected for the reason that SN < SE, V and VN ≃ VE, V. In matrix form MSBAS-2D can be written as

where the matrix consists of the matrix A (described above) multiplied by the east and vertical components of the LOS unit vector S.

The Tikhonov regularization matrix L multiplied by the regularization parameter λ is used for regularizing time series (Tikhonov and Arsenin, 1977). The effect of regularization is similar to applying a low-pass filter and it is recommended when the acquisition times of ascending and descending datasets do not match. The matrix A, as well as the regularization matrices L of zeroth, first and second orders were explicitly derived in (Samsonov, 2010; Samsonov and d’Oreye, 2017). As in the previous case, the unknown east and vertical components of deformation velocities VE and VV for each acquisition epoch and pixel are solved by applying SVD and the 2D surface deformation time series are reconstructed from computed deformation velocities by the numerical integration.

MSBAS-3D

In certain cases, the contribution of the north-south component of displacement is too large to be neglected (SN < SE, V, VN ≃ VE, and VN > VV). By assuming that the motion occurs parallel to the surface, we can however reduce the degrees of freedom from three to two, allowing resolving all three components of motion from (at least) one ascending and one descending DInSAR datasets. This assumption is for instance suited for characterizing gravity-driven processes such as landslides and glacier flows. The constraint can be written in terms of deformation velocities in the following form

where H is the topographic height and and are first derivatives in the north XN and east XE directions. The and values can be computed from DEM using, for example, the Generic Mapping Tools (GMT) grdmath -M DDY/DDX command. Prior low-pass filtering of DEM (e.g. using GMT grdfilter command) is recommended for reducing sensitivity to surface features that do not affect the ground motion originated at some (shallow) depth. It is likely that sliding surface at depth is smoother than DEM. The magnitude of a low-pass filter needs to be sufficiently strong so the resulting surface has a slope equal to an overall slope of the form that host the landslide.

The MSBAS-3D technique that computes constrained 3D surface deformation time series can be written in matrix form as

where the matrix consists of the matrix A (described above) multiplied by the north, east and vertical components of the LOS unit vector S. The matrix represents the motion parallel to the surface constraint. The unknown north, east and vertical components of deformation velocities VN, VE and VV for each acquisition epoch and pixel are solved by applying SVD and the surface deformation time series are reconstructed from the computed deformation velocities by the numerical integration.

For some combinations of acquisition parameters and topographic gradients, (5) may become ill-conditioned, which can be assessed by computing condition number for all possible topographic gradients for given acquisition parameters.

Computational aspects

Matrices in Eqs. (1), (3), (5) typically have hundreds of columns and thousands of rows, corresponding to the total number of SAR images and DInSAR interferograms, respectively. Computation of the Singular Value Decomposition (SVD) of such large matrices is intensive. In case of the 1D analysis, viewing geometry S is irrelevant so SVD computation is performed only once. In case of the 2D analysis, constant viewing geometry S is assumed. This allows performing SVD computation only once, which results in a slight reduction of precision. In case of the 3D deformation analysis, SVD inversion is performed for each coherent pixel in order to account for the spatially varying constraint so varying viewing geometry is implemented as well. The typical processing time of MSBAS-1D and MSBAS-2D is only a few tens of minutes, while the typical processing time of MSBAS-3D can reach several hours (for a typical grid of a few thousand pixels in each dimension) due to the large number of SVD computations required.

Unlike other advanced methods such as StaMPS (Hooper, 2008), MSBAS processing is incremental. The addition of new SAR images requires computing of only the few interferograms satisfying user-defined baseline criteria before executing MSBAS processing.

Test case and data

The ability to unravel surface deformation pattern using the different level of processing of the MSBAS method is illustrated on a large landslide, for which long DInSAR and dGNSS time series, as well as ground truth data, are available (Nobile et al., 2018).

Funu landslide

Funu landslide (2∘31′S, 28∘50′E, Fig. 1) is located in the densely populated City of Bukavu, situated in the Kivu Rift (western branch of the East African Rift System) in eastern Democratic Republic of the Congo. This region is characterised by a high susceptibility to landslides due to a combination of natural predisposing and triggering factors (e.g., high rainfall intensities, moderate to high seismicity, deep weathering and steep landscapes) (Jacobs et al., 2016; Monsieurs et al., 2018; Depicker et al., 2019; Dille et al., 2019). Located in the tropics, the region is affected by persistent cloud cover over most of the year, making DInSAR a key remote sensing tool for the study of landslides in this area (Nobile et al., 2018).

Fig. 1. Funu landslide is located in the city of Bukavu, eastern DR Congo. It lies within the South Kivu Volcanic Province that is formed of sequences of late Miocene to Pleistocene sub-horizontal basaltic lava layers. Outlined in white is shown the area processed with MSBAS. The footprints of the ascending and descending COSMO-SkyMed SAR images used in this study are outlined in black.

This 1.5 km2 slow-moving slope instability was formed thousands of years ago (Moeyersons et al., 2004). The depth of this bedrock landslide can be estimated to be of at least 100 m, i.e. the average height of its main scarp (Dewitte and Demoulin, 2005; Larsen et al., 2010). The landslide developed in sequences of late Miocene to Pleistocene sub-horizontal basaltic lava layers (Fig. 2b) that make up the South Kivu Volcanic Province (Fig. 1). An accurate delineation of the layers does not exist. However, the presence of springs at different elevations shows that their permeability varies, attesting to perched water tables. Weathering of some layers can be very intense with a regolith thickness that can be more than 10 m (Fig. 2c). The distribution of the weathering is highly uneven, sometimes having important changes over short distances of a few tens of meters (Moeyersons et al., 2004). The landslide is currently prone to surface deformation as attested to by continuous damage to infrastructure and buildings (Moeyersons et al., 2004; Balegamire et al., 2017). This landslide is also densely populated, with approximately 80 thousand people living directly on it (Michellier et al., 2016).

Fig. 2. Drone photographs of Funu landslide acquired in September 2018. a) General view of the landslide, with its extent outlined in yellow and green arrows indicating principal direction of movement. The fastest moving region, where discrepancies between DInSAR and dGNSS are observed, is outlined in red. Note constraints from topography and apparent divergent motion at landslides bottom. The total landslide height is about 400 m while the height of the main scarp is up to 100 m. b) Sequences of late Miocene to Pleistocene sub-horizontal basaltic lava layers apparent over the main scarp. c) Gullies delineating the toe of the fastest landslide unit. Their presence highlights the deep weathering of some basaltic layers, with a regolith thickness that can easily exceed 10 m. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Nobile et al. (2018) provided a first analysis of the Funu landslides kinematics and mechanisms by combining one year of DInSAR measurements from 70 Cosmo-SkyMed SAR images acquired between 21 March 2015 and 19 April 2016, processed using the StaMPS SBAS approach (Lanari et al., 2007; Hooper, 2008) and on-site dGNSS measurements. They measured deformation of up to 70 mm/yr in the ascending and descending LOS directions. These measurements were in good agreement with dGNSS, except for one zone of the landslide that corresponds to a recent, smaller subsequent landslide. There, dGNSS measurements showed deformation an order of magnitude larger than DInSAR (up to 500 mm/yr while converted to LOS). The limit in the maximum detectable deformation resolvable by the X-band DInSAR was reached (Nobile et al., 2018). The study by Nobile et al. (2018) confirmed that the Funu landslide is mostly driven by a rotational mechanism with a prominent main scarp and a characteristic back-tilted bench. Several minor scarps, probably linked to secondary shear surfaces, are present in the zone of depletion and delimit individual blocks. From rotational in the depletion zone, the main movement develops a more planar behaviour with flow-like features in its runout and accumulation zones (as suggested from its morphometry). They concluded that the dense urban fabric, the relatively low velocity of the ground deformation, and the general eastward topography with the average slope of about 6.5 degrees (vs 45 degree SAR incidence angle) that avoids shadowing and foreshortening effects make the Funu landslide a particularly interesting site for DInSAR analysis.

COSMO-SkyMed images

We used 176 ascending and 166 descending COSMO-SkyMed images that were acquired between 24 March 2015 and 3 January 2019 (Table 1). Ascending images were acquired with the incidence angle of 41 degrees, range pixel spacing of 1.3 m and azimuth pixel spacing of 2.3 m. Descending images were acquired with the incidence angle of 50 degrees, range pixel spacing of 1.7 m and azimuth pixel spacing of 2.1 m.

Table 1. COSMO-SkyMed SAR data used in this study, θ is azimuth angle and ϕ is incidence angle, N is number of SLC images and M is number of DInSAR interferograms computed for each data set. 2D and 3D time series span 24 March 2015–3 January 2019 interval when both ascending and descending data are available.

Each COSMO-SkyMed set was processed independently with the GAMMA software (Wegmuller and Werner, 1997). Ascending image acquired on 5 February 2017 and descending image acquired on 5 May 2017 were selected as masters for respective sets and remaining images were re-sampled into the master geometry. Multilooking factor of 2 in range and azimuth were used, as well as a maximal perpendicular baselines of 600 m (Fig. 3). The topographic phase was removed using the 1 arc-second (~30 m) resolution Shuttle Radar Topography Mission Digital Elevation Model (SRTM DEM, Fig. 1). Differential interferograms were filtered using the adaptive filter with a filtering function based on local fringe spectrum (Goldstein and Werner, 1998) and unwrapped using the minimum cost flow algorithm (Costantini, 1998). The residual orbital ramp was corrected by applying baseline refinement algorithm implemented in the GAMMA software that utilizes unwrapped phase and topographic height at selected coherent pixels to solve for a refined baseline. Minor interpolation of each interferogram was performed in order to improve the spatial coverage reduced by decorrelation. Then, ascending and descending interferograms were geocoded and resampled using the Geospatial Data Abstraction Library (GDAL) program to a common lat/lon grid with a uniform spatial sampling of 5 m.

Fig. 3. Spatial and temporal baselines of COSMO-SkyMed ascending and descending interferograms used in this study. Brown points mark ascending and descending master images. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Geocoded to a common grid, ascending and descending interferograms with an average coherence after filtering above 0.5 were processed with the MSBAS techniques to produce deformation time series and annual linear deformation rates over Funu landslide and its close surroundings.

dGNSS measurements

The DInSAR-derived deformation rates are validated through a comparison with the dGNSS measurements acquired by repeated measurement campaigns. The dGNSS measurements were captured in a static acquisition mode with Leica GX1230 dual-frequency receivers and post-processed with Leica Geosystems post-processing software. The dGNSS measurements were acquired over already in place benchmarks that are solid concrete structures, such as a water fountains, concrete sewers or street curbs. Measurement duration (5 min) at roving sites was constrained by the safety of operators and equipment in this area, although position accuracy would have benefited from longer sessions. The accuracy achieved for the dGNSS locations was estimated as 7.5 mm for the horizontal and 16.5 mm for the vertical components. These uncertainties refer to the internal precision of the measurement, field conditions, and data post-processing. It was also assumed that all benchmarks and the reference station move at the same rift opening rate of about 2 mm/year (Delvaux et al., 2017).

DEM derivatives

The and values (Fig. 4) were computed from SRTM DEM using the Generic Mapping Tools (GMT) grdmath -M DDY/DDX command. Prior low-pass filtering of DEM was applied using GMT grdfilter command to computed average values in 65 m window and therefore removing surface features that are not expected to affect the ground motion at depth.

Fig. 4. Central first derivatives of 1 arc-second (~30 m) resolution SRTM DEM in north direction (a) and in east direction (b) used in the computation of the three-dimensional deformation with MSBAS-3D.

Results

Fig. 5, Fig. 6, Fig. 7 show annual linear deformation rates for the city of Bukavu computed with the 1D, 2D, and 3D MSBAS techniques. For all methods, the deformation signal of the large Funu landslide clearly stands out.

Fig. 5. Maps of linear deformation rate in the radar line-of-sight (LOS) geometry computed using MSBAS-1D from COSMO-SkyMed ascending (a) and descending (b) SAR data. Displacements along LOS direction have to be interpreted as a one-dimensional component of the full 3-D displacement vector. Red colour represents motion towards satellite and blue colour represents motion away from satellite. 100 m elevation contour lines are plotted in brown. R is reference region. Streets are plotted in gray (downloaded from http://overpass-turbo.eu/). Average SAR backscatter intensity is exposed in decorrelated areas. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
Fig. 6. East (a) and vertical (b) components of linear deformation rate computed using MSBAS-2D from combined COSMO-SkyMed ascending and descending SAR datasets. 100 m elevation contour lines are plotted in brown. R is reference region. Subsiding area that does not deform according to (4) is outlined in red. Streets are plotted in gray (downloaded from http://overpass-turbo.eu/). Average SAR backscatter intensity is exposed in decorrelated areas. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
Fig. 7. North (a), east (b), and vertical (c) components of linear deformation rate computed using MSBAS-3D from combined COSMO-SkyMed ascending and descending SAR datasets. 100 m elevation contour lines are plotted in brown. R is reference region. Subsiding area that does not deform according to (4) is outlined in red. Streets are plotted in gray (downloaded from http://overpass-turbo.eu/). Average SAR backscatter intensity is exposed in decorrelated areas. Points F14, F15, F17, F18, F20 mark location of campaign GNSS sites. For points P1-P6 that approximately correspond to regions of fastest motion time series are provided in Fig. 9. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

The 1D ascending and descending linear deformation rates computed with (1) show distance to satellite (line-of-sight) displacements. Absolute maximal rates of displacement are similar (0.05 m/year) between ascending (Fig. 5(a)) and descending (Fig. 5(b)) datasets, but the extent of the deformation measured over the Funu landslide in the descending map is significantly lesser than in the ascending map.

The 2D east and vertical linear deformation rates computed with (3) show eastward motion with a maximum magnitude of 0.05 m/year (Fig. 6(a)) and downward motion with a magnitude up to 0.02 m/year (Fig. 6(b)). Only the signal of the landslide is shown in the eastward component deformation map, while the vertical displacement map highlights also other deforming regions within the city.

The 3D north, east and vertical linear deformation rates computed with (5) are shown in Fig. 7. The north deformation rate map (Fig. 7(a)) highlights areas with both northern and southern motion with a maximum rate of 0.02 m/year over different places within the city. The east map (Fig. 7(b)) shows eastward deformation with a maximum rate of 0.05 m/year that is nearly identical to the 2D case (Fig. 6(a)) and is dominated by the signal of the large Funu landslide. The vertical map (Fig. 7(c)) highlights only downward motion with a maximum rate of 0.02 m/year. Similarly as with the 2D case (Fig. 6(a)), signals of other deforming regions are captured within the city. The overall scale of the deformation measured with the 3D case is however of lower magnitude.

Annual horizontal displacements measured with MSBAS-3D method are plotted as vectors in Fig. 8, along with the dGNSS deformation vectors. It clearly highlights an advancing landslide lobal front, as well effect of lateral constraints from the surrounding non-moving topography. The directions of movement measured by dGNSS and DInSAR are in general in good agreement at F14, F15, F17, F18 and F20 dGNSS sites, for which the coefficient of determination R2 is greater than 0.75 for north and east deformation components.

Fig. 8. Horizontal deformation rate vectors computed from COSMO-SkyMed ascending and descending DInSAR data with MSBAS-3D technique are in green. Horizontal deformation vectors with 2σ (i.e 95%) confidence interval elipses computed from campaign GNSS data for selected sites F14, F15, F17, F18, F20 are in blue. Funu landslide outline is in red. 100 m elevation contour lines are in white. Background is SAR intensity. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

The 3D surface deformation time series for points P1-P6 (marked in Fig. 7) that approximately correspond to regions of fastest motion within the landslide are shown in Fig. 9. Standard deviation, computed in 5×5 pixels (~25×25 m) window, measures deviation from a mean value rather than precision of the measurement technique because deformation signal and measurement noise cannot be separated. Standard deviation over a larger 33×33 pixels (~165×165 m) window at a region that does not experience observable ground deformation allows measuring precision of this measurement technique. The estimated precision is better than 2.8, 1.4, and 0.5 mm for north, east and vertical components respectively.

Fig. 9. Three-dimensional time series of ground deformation for points P1-P6 (marked in Fig. 7) that approximately correspond to regions of fastest motion. Standard deviation is computed in 5×5 pixels (~25×25 m) window.

Deformation produced by the Funu landslide clearly stands out in all deformation maps. Discrimination between vertical and horizontal movements provided by MSBAS-2D and MSBAS-3D (Fig. 6, Fig. 7) show that, overall, Funu landslide movement is mostly horizontal (to the east). Half-moon deformation patterns highlight rotational movements at the head of the landslide that could not be directly observed from the LOS measures obtained with MSBAS-1D (Fig. 5).

Fig. 5, Fig. 6, Fig. 7 also highlight other ground deformation processes outside of the Funu landslide. At point P2 is observed the deformation associated with another old, deep-seated landslide, as already stressed in Nobile et al., 2018. Deformation observed at the northern limit of the city is probably associated to subsidence of recently urbanized river delta deposits. As similar origin is inferred for documented subsidence over the flat valley bottom east of the study area (red insert). These patterns of deformation were also observed in the field (Kalikone Buzera et al., 2017). Even though the purpose of our research was not to focus on those areas, this further confirms the reliability of the DInSAR processing.

Discussion and conclusions

Working on ascending and descending interferograms computed prior to the MSBAS analysis and geocoded to a common grid, the MSBAS method computed 1D, 2D, and constrained to the parallel surface flow 3D surface deformation time series. From those time series, MSBAS computed map of linear deformation rates for a large, slow-moving, deep-seated landslide in the city of Bukavu. The advantages and disadvantages of each MSBAS and similar methods are summarized in Table 2.

Table 2. Comparison of methods for resolving deformation components from DInSAR and other SAR-derived data.

The LOS measures of deformation obtained with the MSBAS-1D method (Fig. 5) are consistent with the multi-temporal DInSAR SBAS results obtained over a one-year period by Nobile et al., 2018. Overall, the MSBAS results are in agreement with dGNSS measurements with the exception of one zone of the landslide that corresponds to a recent subsequent landslide within the large landslide (Fig. 2b); here, dGNSS measurements show deformation of up to 500 mm/yr (Nobile et al., 2018). These high deformation rates as well as the sharp velocity contrast with the zone directly surrounding it (with a deformation rate nearly an order of magnitude lower), lead to either a loss of coherence or an underestimation of the deformation during the unwrapping of the interferograms. Those are however associated with the intrinsic limitations of DInSAR technique with respect to measure of rapid displacements. Speckle offset analysis performed with the GAMMA software on selected long-timespan pairs did not detect these large deformation rates neither.

The benefits provided by decomposing the deformation components were illustrated in Fig. 5, Fig. 6, Fig. 7, Fig. 8. Compared to MSBAS-2D and 3D, the limitations associated to the traditional 1D LOS analysis (Fig. 5) are clear. While difficult to interpret and communicate to stakeholders unfamiliar with the concept of 1D LOS viewing geometry (Fuhrmann and Garthwaite, 2019), deformation measured on those maps may also lead to misconception about the underlying processes. For instance, individual north, east and vertical deformation components may cancel each other out when the contribution of horizontal components is equal in magnitude but opposite in sign to the contribution of the vertical component. In extreme cases, deformation may appear entirely absent in one (i.e. either ascending or descending) of LOS deformation maps (e.g. Samsonov et al., 2013, Samsonov et al., 2014). The 2D and 3D results (Fig. 6, Fig. 7, Fig. 8) are much easier to interpret thanks to the direct discrimination between horizontal and vertical deformation components. The information gathered about the north component provides additional, non-negligible keys for the understanding of processes at work (see Fig. 8). The 3D solution highlights the presence of an active, advancing landslide front, as well effect of lateral constraints from the surrounding non-moving topography. If this behaviour was anticipated from morphometric analysis of high resolution DEMs and the presence of a large bulge at the landslide toe, this pattern of displacement was completely missed by traditional DInSAR measures (Nobile et al. 2018). With the increase of the understanding of the mechanisms of the landslide (interaction with surrounding topography, future evolution and interaction with downslope stream, etc.), such knowledge also bring information on environmental controls, temporally and spatially-variable external forcing conditions and characteristics of future deformation patterns.

In our example, while the magnitude of the horizontal deformation remains similar between MSBAS-2D and MSBAS-3D, the information about the north component provides clear insights for the understanding of the landslide deformation pattern. The decision of whether to use one or the other will depend on the validity of the Surface Parallel Flow constraint (4) and is case-specific. For many gravity-driven processes, such as landslides and glacier flows, making such an assumption of motion parallel to the surface is a valid approximation. However, even when the assumption of motion parallel to the surface is not strictly correct, the approximate 3D analysis may still provide valuable clues for the understanding of the processes at work that may otherwise go undetected. With our example applied on the city of Bukavu, we also show measure of displacement were expected movement clearly not comply with the assumption of motion parallel to the surface. Besides Funu landslide, the study area encompasses for instance a flat valley bottom that corresponds to an old marsh, recently urbanized and that is subjected to subsidence (Nobile et al., 2018). There, MSBAS-3D incorrectly classifies ground movement as horizontal (Fig. 7(a)), while it is perfectly well grasped with MSBAS-2D (Fig. 6(b)).

In comparison to the method used in Notti et al. (2014); Kalia (2018); Hu et al., 2018, Hu et al., 2019 (Table 2), where the LOS velocity is projected along the steepest downslope to obtain estimates of the landslide motion, MSBAS-3D computes a more general solution. It does not require the motion to occur along the steepest downslope, but only parallel to the surface; which allows to capture lateral motion. The drawback of MSBAS-3D approach is that slope-normal displacements are not allowed. The slope-normal displacements can be produced by mass depletion/accumulation, hydrological loading/unloading, and soil swell/shrink. The magnitude of slope-normal displacements of slow-moving landslides, however, is usually smaller than the magnitude of parallel to the surface displacements. In case of large slope-normal displacements, the Hu et al. (2018) methodology may provide more accurate results.

Other methods commonly used to measure azimuthal deformation have other drawbacks (Table 2). Speckle tracking (e.g. Wright et al., 2004; De Michele et al., 2010; Raucoules et al., 2013; Shi et al., 2018) or Multiple Aperture Interferometry (MAI) (e.g. Gourmelen et al., 2011) approaches are for instance best suited for measuring meter-scale deformation (Fuhrmann and Garthwaite, 2019). When the measured deformation is of lower magnitude (e.g., millimetre- to decimetre-scale displacement associated to slow moving landslides (Hungr et al., 2014), the north-south displacement component accuracy provided is much lower to what can be reached with conventional DInSAR. Also, the competing techniques presented in the introduction commonly work on individual (or a series of disconnected) interferograms to derive the 3D deformation. By potentially missing controls of temporally and spatially-variable external forcing on the process kinematic, they are more adapted to the study of sudden deformation events such as earthquakes or volcanic eruptions (Hu et al., 2014). Here lie the strengths of the MSBAS method: producing 3D deformation time series with the accuracy of conventional DInSAR and the spatio-temporal resolution of modern SAR sensors, it can provide unrivalled insights for studying both long and short-term kinematics of surface processes. Those benefits associated to the 3D decomposition of surface deformation were demonstrated here on a large slow-moving landslide, but this method can be applied to studying many other gravity-driven processes producing motion parallel to the surface.

Authors: Sergey Samsonov, Antoine Dille, Olivier Dewitte, François Kervyn, Nicolas d’Oreye

admin

admin

One thought on “Satellite interferometry for mapping surface deformation time series in one, two and three dimensions: A new method illustrated on a slow-moving landslide

Leave a Reply

Your email address will not be published. Required fields are marked *