Monitoring activity at the Daguangbao mega-landslide (China) using Sentinel-1 TOPS time series interferometry

The Daguangbao mega-landslide (China), induced by the 2008 Wenchuan earthquake (M w = 7.9), with an area ofapproximately 8km 2 ,isoneofthelargestlandslidesintheworld.Expertspredictedthatthepotentialriskand instabilityofthelandslidemightremainformanydecades,orevenlonger.Monitoringtheactivityofsuch alarge landslide is hence critical. Terrain Observation by Progressive Scans (TOPS) mode from the Sentinel-1 satellite provides us with up-to-date high-quality Synthetic Aperture Radar (SAR) images over a wide ground coverage (250×250km),enablingfullexploitationofvariousInSARapplications.However,theTOPSmodeintroducesaz-imuth-dependent Doppler variations to radarsignals, whichrequiresanadditional processing step especially for SARinterferometry.Sentinel-1TOPSdatahavebeenwidelyappliedtoearthquakes,buttheperformanceofTOPS data-based time series analysis requires further exploitation. In this study, Sentinel-1 TOPS data were employed to investigate landslide post-seismic activities for the ﬁ rst time. To deal with the azimuth-dependent Doppler variations,aprocessingchainofTOPStimeseriesinterferometryapproachwasdeveloped.SincetheDaguangbao landslideisasaresultofthecollapseofawholemountaincausedbythe2008Mw7.9Wenchuanearthquake,the existingDigitalElevationModels(DEMs,e.g.SRTMandASTER)exhibitheightdifferencesofupto approximately 500 m. Tandem-X images acquired after the earthquake were used to generate a high resolution post-seismic DEM. The high gradient topographic errors of the SRTM DEM (i.e. the differences between the pre-seismic SRTMandtheactualpost-seismicelevation),togetherwithlowcoherenceinmountainousareasmakeitdif ﬁ cult to derive a precise DEM using the traditional InSAR processing procedure. A re-ﬂ attening iterative method was hence developed to generate a precise TanDEM-X DEM in this study. The volume of the coseismic Daguangbao landslide was estimated to be of 1.189 ± 0.110 × 10 9 m 3 by comparing the postseismic Tandem-X DEM with the preseismic SRTMDEM,which is consistent with the engineering geological survey result. The time-series re-sults from Sentinel-1 show that some sectors of the Daguangbao landslide are still active (and displaying four sliding zones) and exhibiting a maximum displacement rate of 8 cm/year, even eight years after the Wenchuan earthquake. The good performance of TOPS in this time series analysis indicates that up-to-date high-quality TOPS data with spatiotemporal baselines offer signi ﬁ cant potential in terms of future InSAR applications.


Introduction
At 6:28 UTC On 12 May 2008, a devastating Mw 7.9 earthquake hit the Sichuan province in southwest China,resulting in 69,227 deaths,18,195 missing,374,216 injured and 5,362,500 buildings collapsed (Cui et al., 2009).In addition to the immediate devastation caused by shaking, the earthquake triggered more than 60,000 destructive landslides over an area of 35,000 km 2 (Huang and Fan, 2013) that created more than 30 landslide dams that threatened both downstream and upstream areas (Chigira et al., 2010).Millions of cubic meters of loose soil, debris and rocks were accumulated on the hill-slopes within the area of the Wenchuan earthquake.This debris subsequently became a source of material with the potential to lose stability during wet seasons due to rainfall infiltration, and thus quickly evolve into deadly landslides (Zhang et al., 2014).Strong earthquakes (e.g.Mw 7.0 or higher) not only trigger coseismic landslides, but can also provide loose material or lead to increased post-seismic slope instability over extended periods of time (e.g. 1 to 10 years) (Tang et al., 2011).How long it will take for the landslide/debris flow frequency to return to pre-earthquake levels depends on a range of factors, including rainfall intensity, natural re-vegetation, self-stabilization processes on slopes, erosion, and valley incision (Huang and Fan, 2013).Huang and Fan (2013) anticipated that debris flows directly resulting from sediment movements during the 2008 earthquake may remain active for another two decades.Therefore, it is crucial to monitor the activity of post-seismic landslide/debris flow hazards in this region.
The Daguangbao landslide, induced by the 2008 Wenchuan earthquake with a volume of 1.167 × 10 9 m 3 (Huang and Fan, 2013), is one of the largest active landslides in the world (Table 1).This immense unprecedented landslide has drawn much attention from geologists and geophysicists.Huang et al. (2011) investigated the characteristics, failure mechanism and dynamic processes of the landslide.Xu et al. (2013) analyzed the characteristics of the formation mechanism and kinematics of the landslide.Chen et al. (2014) quantified the mass wasting volume of the landslide using topographic change information derived from InSAR data.However, all of these previous studies only focused on characterizing and modeling the co-seismic period of the Daguangbao landslide and, as such, nothing is known about its current activity (e.g.whether a new equilibrium could be initiated by continuous or intermittent movements).Furthermore, monitoring small displacements taking place over the entire area of the post-seismic landslide would provide valuable information on the kinematics of the loose material accumulated by the Wenchuan earthquake, but would be a time-consuming and difficult task.
Synthetic Aperture Radar Interferometry (InSAR) technology provides the capability to monitor precise surface displacements over time with a wide coverage in a time-efficient and labor-saving manner (Tomás et al., 2014a).Time Series InSAR (TS-InSAR) technology, characterized by its reduction of de-correlation and atmospheric effects, has been used successfully in various applications and fields of research as city subsidence (Chen et al., 2016;Dai et al., 2015;Herrera et al., 2010b;Jiang et al., 2011;Kim et al., 2015;Zhang et al., 2012), earthquakes (Wright et al., 2004;Feng et al., 2015;Feng et al., 2016;Wang et al., 2014) and landslides (Chen et al., 2014;Herrera et al., 2010a;Herrera et al., 2013;Singleton et al., 2014;Tomás et al., 2014b).European Space Agency (ESA) Sentinel-1 satellites have provided us with upto-date, high-quality images covering a wide area in the Terrain Observation by Progressive Scans (TOPS) mode.TOPS mode SAR images are comprised of bursts characterized by large Doppler centroid variations, which is significantly different from those in the strip map mode.Residual azimuthal ramps appear during the interferometric phase where azimuth mis-registration errors greater than a few thousandth per sample occur (Lanari et al., 2015).There have been more and more research works and technical reports on Sentinel-1 TOPS data interferometry (Prats-Iraola et al., 2015;Yague-Martinez et al., 2016) with its main application to earthquakes (Feng et al., 2016;Grandin et al., 2016;Wen et al., 2016).The performance of the TOPS data-based time series interferometry and its application to landslides thus requires further exploitation.
In this paper, we use the TS-InSAR integrated Atmospheric Estimation Model (TSInSAR_AEM) method to monitor the current activity of the Daguangbao landslide with TOPS Sentinel-1A images.TanDEM-X images are employed to generate a high-resolution postseismic Digital Elevation Model (DEM), which is then used to simplify interferometirc processing and time series analysis.One 3 m multi-spectral Dove image (with R, G, B and NIR bands) is utilized to asssist with our interpretation and analysis of the landslide.

The geological setting and evolution of the Daguangbao landslide
The study area is located in the Tibetan and Qiang Autonomous Prefecture, southwest China (Fig. 1a).Fig. 1b and c show the geological setting of the Daguangbao landslide subsequent to the coseismic movement of the 2008 Wenchuan earthquake.The landslide is located at the north flank of the Dashuizha anticline, in a mountainous zone with elevations between 1000 and 4500 m a.s.l. and is mainly NE oriented.The lithology of the outcrop in the study area is mainly carbonated rocks from the Precambrian to Cretaceous age (Huang et al., 2011;Huang et al., 2014;Zhang et al., 2015) and comprises: (a) the Feixianguan group of the Triassic system (Tf) which consists of violet siltstone, silty mudstone with a small amount of shell and micritic limestone, (b) the Liangshan group (Pl), Yangxin group (Py1), Longtan group (Py2) and Wujiaping group (Pw); these groups belong to the Permian system which is composed of medium to thick layered crystalline and micritic limestone with flint, (c) the Zongchang group of the Carboniferous system (Cz) which consists of violet red sandy mudstone, (d) the Shawozi group of the Devonian system (Ds and Ds) which contains dolomite with phosphate mineral rocks, (e) and the Dengyin group of the Sinian system (Zd and Zg) which is composed of limestone, interbeds of red mudstone and dolomitic rocks.
The 2008 Wenchuan earthquake generated an intensive surface motion causing maximum peak ground accelerations of 957.7 cm/s 2 (EW), 652.9 cm/s 2 (NS) and 948.1 cm/s 2 (UD) at the Wolong station located 23 km from the epicenter, and endured for a total of 120 s (Chen et al., 2014), and triggered the Daguangbao landslide, with a maximum coseismic displacement of about 1.5 km (Huang and Fan, 2013).The Daguangbao landslide, with an estimated affected area of 7.3-10 million m 2 and a volume of 1.167 × 10 9 m 3 , was the largest landslide induced by this earthquake, originated in the Longmen Shan fault zone at the eastern margin of the Tibetan Plateau (Huang and Fan, 2013).It remains unknown whether the landslide has become stable after the seismic event or whether it will remain unstable until a new equilibrium is achieved, which will be addressed in this paper.
The landslide has an open V-shape sliding surface in a plan view of more than 1 km in length (Fig. 2a), defined by two main discontinuities corresponding to the bedding (dipping towards the N 30-35°) and a fault that trends SW-NE and dips 75°SE (Chigira et al., 2010).The source area (i.e. the depletion area according to Cruden and Varnes (1996)) is approximately 2.4 km long and 1.2 km wide (Huang et al., 2011) (Fig. 1).The accumulated mass has formed a hummock deposit of 3.2 km in length and 2.2 km in width (Huang et al., 2011) with a maximum thickness of 521 m (Huang et al., 2014) (Fig. 2); this deposit is mainly composed of loose materials generated by the collapse of the slope and the subsequent sliding processes.Furthermore, a large fractured non disaggregated rock blocks of almost 1 km in length (and which partially preserves its original structure) is located at the toe of the landslide.These loose deposits may be in a quasi-unstable situation.Fig. 2c and d show two photos of the Daguangbao landslide taken on 19 October 2012, whose points of view are labeled in Fig. 2a.It should be noted that the debris deposits consist of bare rock and loose rock boulders, which act as strong scatterers in SAR images and provide high coherence in InSAR measurements.

Satellite data
Sentinel-1A, a space mission of the Europe Space Agency of the Copernicus Programme, was launched on 3 April 2014.This mission provides C-band continuity following the retirement of ERS-2 and the end of the Envisat mission (ESA, 2016).Compared to other SAR missions, its short revisit time (twelve days on its own, and six days with Sentinel-1B) and short baselines (i.e.orbital separations) make Sentinel-1A promising in terms of InSAR applications.Its new type of ScanSAR mode, known as TOPS SAR, has a similar pixel resolution (about 3.5 m in the range direction and 14 m in the azimuth direction) to ScanSAR but with an enhanced signal-to-noise ratio distribution (Farkas et al., 2015).It also has wider coverage (up to 250 km) than previous SAR sensors.This is achieved by a rotation of the antenna in the azimuth direction, and requires very precise coregistration because even errors under pixel accuracy can introduce azimuth phase variations caused by differences in Doppler-centroids (Farkas et al., 2015).Sentinel-1 SAR data have been regularly acquired over China since March 2015.
Ascending (satellite moving north) Sentinel-1A images have a more favouriate geometry than descending (satellite moving south) ones for the Daguangbao landslide, and the number of the available descending images is limited.Hence only 15 ascending scenes of Sentinel-1A TOPS mode Interferometric Wide Swath (IWS) acquisitions between March 2015 and March 2016 were used in this study (Table 2).As evidenced in Fig. 3, The Sentinel-1 orbit maintenance strategy ensures a precise ground-track repeatability and hence small InSAR baselines (i.e.orbital seperations) on the order of 150 m (Yague-Martinez et al., 2016).Therefore, 53 interferograms with high coherence were generated, providing us with sufficient redundant observations to ensure the quality of final time series results.
For the Daguangbao landslide, the existing DEMs (i.e.SRTM and ASTER DEMs) were both acquired prior to 2009 (before the earthquake).As mentioned earlier, the coseismic landslide caused the collapse of a whole mountain.As a result, the differences between the SRTM DEM and the actual postseismic elevation (commonly termed as topographic error in interferometric processing) should represent the landslide thickness, with a maximum thickness of 521 m (Huang et al., 2014).In this study, the Tandem-X (TerraSAR-X add-on for Digital Elevation Measurement) data acquired on 21 July 2011 were used to generate a high resolution DEM covering the Daguangbao landslide.TanDEM-X (the first bistatic SAR mission) is formed by adding a second, almost identical, spacecraft to TerraSAR-X and flying the two satellites in a closely controlled formation at a typical distance between 250 and 500 m (German Aerospace Center, 2010).This unique twin satellite constellation allows the generation of the DEM with unprecedented accuracy, pole-to-pole coverage and uniform quality according to the  HRTI-3 specifications (EoPortal Directory, 2015).In the Tandem-X interferometric processing in this study, phase unwrapping errors caused by high gradient topographic errors (i.e. the differences between SRTM and the actual elevation) and low coherence made it a challenge to derive a precise DEM.Therefore, a re-flattening method was developed to address the high gradient topographic errors; details are presented in Section 4. In order to interpret InSAR time series results, a multi-spectral image (with R, G, B and NIR bands) with a 3 m spatial resolution from Planet Labs, Inc. (PL) was utlised.These up-to-date and high resolution images play an important role in assisting us to better interpret InSAR results and understand the evolution of the Daguangbao landslide.

TOPS data coregistration
Due to the steep azimuth spectrum ramp presented in each burst, TOPS interferometry requires highly stringent coregistration, especially in the azimuth direction.This is because the implementation of TOPS (obtained by steering the antenna from aft to fore within a burst) introduces azimuth-dependent Doppler variations on the received signal (Nannini et al., 2016).The residual azimuthal ramps in the interferometric phase are presented even when azimuth mis-registration errors are greater than a few thousandths of a pixel (Lanari et al., 2015).Therefore, coregistration is a key step towards TOPS interferometric processing.
The Sentinel-1 IWS SLC product can be acquired from the ESA's Sentinel Scientific Data Hub website.This product includes three subswaths ("burst SLC") which are obtained by processing the bursts over it.As shown in Fig. 4, the whole SLC is covered by several bursts and the overlaps appear in individual bursts in both azimuth (between sub-sequent bursts) and range (between neighboring sub-swaths).Following the radiometric calibration, the individual bursts can be merged into one SLC (Single Look Complex).The generated mosaics SLC are typically seamless in both range and azimuth based on the very high geometric and radiometric standards of Sentinel-1 (Wegmüller et al., 2015).Subsequent to this, the equation used to resample one SLC image into the geometry of a reference SLC image (rough coregristration) is derived with a consideration of the terrain topography.Coregistration accuracy can be reached at a 1/100 azimuth pixel level by implementing the cross correlation matching procedure (potentially iteratively performed).The strong Doppler variations within each burst result in significant phase jump effects, even with coregistration accurcy at 1/100 pixel level in the azimuth direction.
A spectral diversity method (Scheiber and Moreira, 2000) that considers the double difference phase in the burst overlap regions is performed (potentially iteratively) for further refinement, in order to cope with phase jump issues between sub-sequent bursts.The Doppler Centroid runs through a relatively steep spectral ramp and the azimuth spectrum deramping process is applied to remove it.Finally, the coregistrated SLC is acquired with 5/1000 azimuth pixel accuracy, and can be regarded as the 'normal' SLC in the subsequent processing steps.In this study, these processes were largely implemented using the GAMMA software.

TSInSAR-AEM time series analysis
InSAR time series has been developed over the past decade, and can reduce various decorrelations by combining multiple SAR images.Persistent Scatterers Interferometry (PS-InSAR) (Ferretti et al., 2001) and Small Baseline Subsets (SBAS) (Berardino et al., 2002) form the two most widely used InSAR time series methods.In this study, a modified SBAS method (TSInSAR-AEM), as proposed by Li et al. (2009), which integrates the Atmospheric Estimation Model, was used to process the available SAR images.A brief description of this method is provided in this section.
Assuming a set of N + 1 co-registered SLC images acquired at the ordered time (t 0 , t 1 , …, t n + 1 ), the first step is to generate and optimize image pair combinations.In this study, this was achieved by the constraints on the small spatial and temporal baseline, as well as the frequency shift between the Doppler centroids.In some cases, the constraints would lead to several seperate subsets of interferograms.Certain specific interferograms should be formed to link these subsets.
For these M small baseline interferograms, convertional interferometric processing was applied and the flattening and topographic effects were removed.Phase unwrapping was performed in the 2D space for each interferogram.In this case, we used a minimum cost flow unwrapping (MCF) algorithm (Chen and Zebker, 2002) to acquire the phase changes in the radar line of sight (LOS).Subsequently, a phase closure technique (Biggs et al., 2007) was used to identify the major unwrapping errors remaining in the interferograms.Subsequent to this, for an arbitrary coherent pixel (x,y) in interferogram i (generated by SAR images acquired at time t 1 and time t 2 ), its value can be modeled as, where φ disp , φ topo , φ atm and φ res denote the phase components along the radar Line of Sight (LOS) due to land deformation, topographic error, atmospheric artifacts and decorrelation/thermal noise, respectively.λ is the wavelength and θ is the incidence angle.d t2 (x,y) and d t1 (x,y) are the cumulative displacements (with respect to reference time t 0 ) along the LOS at time t 2 and t 1 , respectively.△Z(x,y) is the topographic error between the actual elevation and the DEM used in the interferometric processing.This technique is performed by way of a pixel-by-pixel temporal analysis.Accordingly, the dependence of Eq. (1) on the (x,y) variable is not explicitly mentioned hereafter.Therefore, for a given pixel in the time span, Eq. ( 1) can be generalized into a matrix equation as, where, and D is the Atmospheric Phase Screen (APS) variable for the specific date t i .t m and t s denote the master and slave images for one interferogram, while Vector A for t i is given by, Based on the fact that the atmospheric signal phase component is highly correlated in space but poorly in time (Ferretti et al., 2000), a known temporal deformation model which is dependent on different cases can be applied to V and helps to separate atmospheric effects from deformation signals.If no temporal deformation model is known, a temporarily linear velocity (TLV) model can be employed to help estimate theφ atm at time t i , It is clear that the TLV model only makes an assumption on two consecutive intervals and thus has a minimum impact on the whole TS (Li et al., 2009).Eqs. ( 6) can be rewritten as a matrix format, as follows, Combining Eqs. ( 2) and ( 7), the final equation would be, If the spatial baseline network is well connected (i.e.there are no separate subsets), we will have N ≥ S, and its solution can be obtained in a least squares sense.In this study, the algorithm was applied iteratively until convergence was achieved.The details can be found from Li et al. (2009) and Zhou et al. (2016).
4. Tandem-X data processing and results analysis 4.1.Re-flattening iterative method SRTM DEM (1-arc-second) was used in this case, in which the height reference datum was transformed from a geoid datum into an ellipsoidal datum (i.e.WGS 84) as suggested in Li et al. (2013).The original Tandem-X images acquired were in the Coregistered Single-look Slant-range Complex (CoSSC) format.So the basic steps initiate with the two co-registrated data.As with the standard introfereometric processing procedure, the interferogram is generated before the flattening and topographic effects are removed by way of the external SRTM DEM.As a result of using the bistatic mode in Tandem-X, certain stepssuch as interferogram flattening, phase calibration, phase simulation and geocodingare slightly different (Pandit et al., 2014).The adapted filter (i.e. the Goldstein filter; Goldstein and Werner (1998)) is applied to the interferogram for the purpose of noise reduction.The filtered interferogram is further unwrapped by way of the triangulation based minimum cost flow method, and all the pixels with coherence lower than a threshold of 0.2 were masked.Thereafter, the unwrapped interferogram is used to generate a height map using the precise baseline information, which can be refined by external ground control points (GCPs) or an existing DEM.In the final step, the height map is geocoded to the geodetic coordinate system relative to the WGS84 ellipsoid.However, this standard processing did not result in a good DEM in this case, mainly because of the phase unwrapping errors brought about by the high gradients of topographical errors (i.e. the differences between the pre-seismic SRTM and the actual elevation after earthquake) and the low coherence in the mountainous areas.As such, a re-flattening iterative method with two refinement steps (i.e.refining-orbit and re-flattening) was developed in this study (Fig. 5).

Refining-orbit
In this study, GCPs were utilized for the purpose of precise orbit estimation.Without external GCPs data, the common method is to extract the height from the existing DEM at certain uniform distributed points as GCPs.However, in this case there was significant topographic errors for the external DEM, as it was acquired prior to the earthquake.We therefore propose the introduction of a refining-orbit step into the process.After the first phase to height conversion, the height results on the GCPs (based on pixels) were used for comparison purposes with the DEM and only points with a difference smaller than 2 m remained.A non-linear least-squares solution algorithm was used to determine the baseline parameter values at these points.The phase to height conversion was then performed once again, based on the new set of precise baseline parameters.

Re-flattening
As mentioned earlier, the topographic errors (i.e. the differences between the pre-seismic SRTM and the actual elevation after earthquake)  are up to approximately 500 m within a horizontal distance of less than 4 km.Due to the phase unwrapping error caused by the high gradients of topographic errors and the low coherence in the surrounding areas of the Daguangbao landslide caused by the dense vegetation, there is no guarantee for the quality of the resultant DEM in the first round.As shown in Fig. 6b and c, when the derived DEM was used in the next iterative interferometric processing, residual fringes were evident.As the Tandem-X was acquired at the same time (i.e.there should be no deformation), the residual fringes were dominated by the topographic errors, resulting in the imprecision of the resultant DEM.It is clear that the residual fringes gradually decrease with the iteration (from Fig. 6a to d), suggesting an increase in the precision of the DEM.We accepted the results of the third iterative interferometry process, since no topographydependent residual fringe is observed in the 3rd re-flattened interferogram (Fig. 6d).

Differences between tandem-X and SRTM DEMs
The postseismic Tandem-X DEM was generated by way of the abovementioned re-flattening iteritive method.The SRTM4, acquired in 2001 (before the earthquake), was interpolated into the same grid for the purpose of comparison.EoPortal Directory (2015) reported that Tandem-X DEMs have a relative vertical accuracy of 4 m for areas with a slope greater than 20%, and Jet Propulsion Laboratory (2016) suggested that the relative vertical accuracy of SRTM4 is smaller than 10 m.Two recent validation studies suggested that, in mountainous areas in China (e.g. the Three Gorges region), the absolute vertical accuracy of SRTM4 is 15.08 m (Li et al., 2015) whilst the absolute vertical accuracy of Tandem-X DEM is about 4.9 m (Li et al., in preparation), suggesting an uncertainty of 15.8 m in their height differences.These figures should be representative for the study site in this study since it is also in a mountainous area in China.The difference map between Tan-DEM-X DEM and SRTM4 over the study site is illustrated in Fig. 7c with a maximum difference of greater than 500 m and a rough uncertainty of 15.8 m.Two significant variations are evident in the middle section, which demonstrates once again that the Daguangbao co-seismic landslide occurred as a result of the 2008 Wenchuan earthquake.It is clear in Fig. 7 that the collapsed moutain produced elevation changes greater than 400 and 500 m in the zones of depletion and accumulation of the coseismic landslide, respectively.
After we acquired the postseismic Tandem-X DEM, the volume of the coseismic landslide can be estimated by comparing it with the preseismic SRTM DEM.Based on the height differences of these two DEMs, the total volume of the landslide accumulated mass was estimated to be 1.189 ± 0.110 × 10 9 m 3 .It is clear that the volume of estimated landslide accumulated mass is consistent with the value derived from recent studies (Table 3).

TOPS time series results and analysis
With the TOPS datasets and the Tandem-X DEM, the time series results were acquired using the TSInSAR-AEM method, as mentioned previously.The landslide surface of Daguangbao, which is mainly composed of bare rock mass and loose rock boulders (Fig. 2c and d), provides good coherence for InSAR data, and ensures the high quality of our TSInSAR-AEM performance.Moreover, the high number of observations taken from the 53 interferograms also guarantees the precision of the time series results.
Fig. 8 shows the annual displacement velocity of the Daguangbao landslide in the LOS direction.The Tandem-X DEM is shown in the 3D elevation wireframe as the bottom layer.Four main sliding zones exhibit "very slow" velocities (IGUS/WGL, 1995), with a maximum displacement rate of 8 cm/year within the whole Daguangbao landslide; these main sliding masses are labeled as SM-a, SM-b, SM-c and SM-d in Fig. 8.The southeast corner of the landslide is the location of the lowest area of the landslide.As such, this section exhibits positive displacement rates which means that there are debris accumulations.
The relationships between line of sight (LOS) and downslope displacements are summarized in Fig. 9 (note that the displacement is assumed to be in the downslope direction).When the slope faces the sensor, downslope displacements will result in the ground surface moving towards the sensor in the line of sight if its slope angle is smaller than the incidence angle (Scenario 1); the sensor will not be able to capture downslope displacements if the slope angle is equal to the incidence angle (Scenario 2); downslope displacements will lead to the ground surface moving away from the sensor in the LOS direction if the slope angle is greater than the incidence angle (Scenario 3, e.g.SM-b in this study).When the slope faces away from the sensors, downslope displacements will result in the ground surface moving away from the sensor if the slope angle is equal or smaller than (90 −θ) (Scenarios 5 and 6, e.g.SM-a, SM-c, SM-d in this study); the sensor will not illuminate the ground surface (i.e.shadow) if the slope angle is greater than (90 − θ) (Scenario 4).Fig. 10 shows the time series result of the whole landslide during the period from March 2015 to March 2016.The dates of the images and the days relative to the reference date (20150331) are presented, from which the spatial pattern of surface ground deformation consistent with the reactivation of some parts of the deposits accumulated by the coseismic landslide is evident.
The main geomorphological features of the four recognized sliding masses were mapped based on geomorphological photo-interpretation of visible landslide features through a multi-spectral 3 m spatial resolution image, PL (Fig. 11a), which was acquired on 30 August 2015.Only the main scarps can be recognized and mapped as it is difficult to be precisely delineated due to their incipient deformations, the continuous   [4]   Note: SRTM DEM was used as the preseismic DEM, whilst the postseismic DEM was derived from:  west of this zone because of the shadow effect caused by the SAR oblique observation geometry.The available results of SM-a indicate that even eight years after the earthquake, this area, located on the main scarp of the co-seismic landslide in which more than 250 m thickness of loose materials (Fig. 2b) were accumulated during the event, remains active.These displacements are apparently associated with two sliding masses (SM-a1 and SM-a2 in Fig. 11) composed of debris (debris slide, according to Cruden and Varnes, 1996) that present different downslope directions.
Slumping masses SM-b1 and SM-b2 affect the loose deposits accumulated on the east slopes placed at the toe of the coseismic landslide.These deposits with a thickness ranging from 150 to 550 m (Fig. 2b) were accumulated as a result of the coseismic landslide over the opposite Huangdongzi valley wall and thus, their instability is evident.A clear scarp can be recognized on the head of SM-b2 (Fig. 11) and according to the slope aspect, two different debris slides (SM-b1 and SM-b2) with a maximum displacement rate of 7.6 cm/year are recognized.
On the SE sector of the coseismic Daguangbao landslide there exist a main scarp.Behind this scarp there is a secondary scarp that clearly delimits the area of maximum LOS displacements (up to 6.6 cm/year) of the debris slide.The accumulated deposits reach in this area up to 500 m of thickness (Fig. 2b) and the movement of the slide presents a SE component that seems to follow the Huangdongzi valley direction.
A debris slide with a minor extension and lower displacement rates (SM-d in Fig. 11) has been recognized in the northeast part of the coseismic landslide.This area, in which up to 350 m of loose deposits were accumulated as a result of the coseismic landslide, presents a NE movement up to 5.4 cm/year and is associated with an obvious scarp (Fig. 11).
Consequently, it can be concluded that the active landslides identified by means of InSAR consist of independent debris slides (according to Cruden and Varnes (1996)) contained within the main coseismis landslide and are mainly associated with the precarium equilibrium of high thicknesses of loose materials deposited by the landslide induced by the 2008 Wenchuan earthquake in a metastable state.Rainfall is a possible factor that controls the activity of these postseismic landslides and analysis of long time series of InSAR data will allow for the assessment of this assumption in the future.

Conclusions
In this paper, the current activity of the immense Daguangbao landslide has been monitored by way of Sentinel-1 TOPS time series interferometry.In order to solve the azimuth ramp between sub-sequent bursts caused by azimuth-dependent Doppler variations in TOPS data, special attention has been paid to their image coregistration.As the interferometric processing with TOPS data is significantly different from that with strip map data, the whole processing chain for TOPS data coregistration and interferometry has been developed and presented in detail, in which the cross correlation matching and the spectral diversity methods are iteratively used until accuracy of 0.001 pixel in the coregistration (especially in the azimuth direction) is achieved.This interferometric processing chain ensures that the Sentinel-1 interferograms are not affected by the azimuth-dependent Doppler variations.
Moreover, it should be noted that, as a result of the 2008 Wenchuan earthquake, the topographic error of the pre-seismic DEM extended to 500 m, introducing a potential error source to the time series analysis, which is another challenge in this study.In order to resolve this error, we used Tandem-X data to derive a high resolution post-seismic DEM.The high gradient topographic error and the low coherence in the mountainous area made it difficult to acquire accurate results.An iterative re-flattening method was thus developed to improve the generated DEM, in which the refine-orbit and re-flattening process were performed.Using this method, a high resolution Tandem-X derived DEM was acquired.Compared with the previous DEM, the co-seismic motion of the Daguangbao landslide was clearly presented.
In this study, fifty three high quality interferograms were acquired from 15 images of the Daguangbao landslide, which provided sufficient redundant observations and ensured the quality of the final results.The following three factors contributed towards the high quality achieved: (1) high coherence from the debris surface of the landslide; (2) the small perpendicular baselines of Sentinel-1 images (the baselines of almost all the interferometric pairs were smaller than 150 m); and (3) the short revisit time of the Sentinel-1 system.Based on these high-quality interferograms the advanced TSInSAR-AEM method was applied to extract InSAR time series results.From the time series results, four main debris slide zones were identified within the area of the coseismic Daguangbao landslide and the maximum displacement rate was found to extend up to 8 cm/year, indicating that some parts of the mass slid during the Wenchuan earthquake remain active even eight years after the Wenchuan earthquake with a "very slow velocity" due to their metastable state of equilibrium.
It is clear from this study that, with small temporal and spatial baselines, Sentinel-1 A can provide a large number of interferograms of high quality and shows great potential for InSAR time series applications.The recent launch of Sentinel-1B satellite improves the revisit time of Sentinel-1 from 12 to 6 days, which is expected to further promote the use of Sentinel-1 data for land monitoring and emergency services.

Fig. 1 .
Fig. 1.(a) Location of the Daguangbao landslide.(b) Geological settings of the study area.(c) Cross-section of the Daguangbao landslide (adapted from Zhang et al. (2015) and Huang et al. (2014)).The two main discontinuities that define the sliding surface are labeled in the map as bd (bedding) and f1(fault).

Fig. 2 .
Fig. 2. (a) Plan view of the landslide.(b) Loose deposits thickness located in the Huangdongzi valley after the Wenchuan earthquake (Huang et al., 2014).(c) Main scarp and loose deposits.(d) Downslope general view of the landslide.Note that (i) the two main discontinuities that define the sliding surface are labeled in the map as bd (bedding) and f1 (fault) in Fig. 1a, and (ii) the displaced mass preserving its original structure and the lower part of the south flank of the landslide in the center and on the right side of the picture, respectively.D: loose material deposits; TC: talus cones deposits; MS: main scarp; Sfl: south flank; Nfl: north flank.

Fig. 3 .
Fig. 3.The spatial and temporal baselines of the Sentinel-1 data used in this study.

Fig. 6 .
Fig. 6.Interferograms (after removing flat and topographic effects) for each iteration.The white lines indicate the boundary of the coseismic landslide (in SAR Range-Doppler Coordinate System).

Fig. 5 .
Fig. 5. Flow chart for the 're-flattening' iterative method used for DEM generation.

Fig. 7 .
Fig. 7. Differences between SRTM and Tandem-X DEMs.Note that in (a) and (b) the black dashed lines indicate the boundary of the coseismic landslide (in SAR Range-Doppler Coordinate), and in (c) positive and negative values indicate mass loss and mass increase, respectively.Note that grey in (b) indicates areas with invalid values due to the shadow of SAR geometry or low coherence caused by dense vegetation.

Fig. 8 .
Fig. 8.The LOS displacement rate map of the Daguangbao landslide obtained from the Sentinel-1A data (negative values indicate the ground surface is moving away from the sensor).Black arrows indicate the maximum topographical gradient (steepest slope) directions derived from the postseismic Tandem-X DEM.
Fig. 10.Time series results (from 24/04/2015 to 25/03/2016) of the whole landslide.The colour bar is identical to that in Fig. 8 and the magenta line indicate the boundary of the Daguangbao landslide.

Fig. 9 .
Fig. 9. Relationships between the line of sight (LOS) and the downslope displacements for different slope orientations (Scenarios 1, 2 and 3: facing the sensor, and Scenarios 4, 5 and 6: facing away from the sensor) and slopes (α and α′).

Fig. 11 .
Fig. 11.(a) High resolution PL Dove true-colour imaging of the Daguangbao landslide.(b) InSAR-derived LOS displacement rate map (the negative value indicates the target moved away from the sensor).

Table 2
Detailed parameters of the Sentinel-1 data.

Table 3
Comparisons of the landslide accumlated mass volumes.