Recent Landslide Movement in Tsaoling , Taiwan Tracked by TerraSAR-X / TanDEM-X DEM Time Series

The Tsaoling Landslide in Taiwan has captured attentions of researchers worldwide due to its frequent catastrophic failure and distinctive features. Thanks to the launch of TerraSAR-X/TanDEM-X (TSX/TDX) constellation, retrieval of global DEM with high spatial resolution and accuracy becomes possible, which is extremely useful for the study of natural hazards (e.g., landslides) globally. We attempt here for the first time to track recent landslide movements in Tsaoling Taiwan by analyzing DEM time series reconstructed from TSX/TDX image pairs. Quality improvement of InSAR derived DEM through an iterated differential operation is addressed during the data processing. Five cliffs and the Chingshui River are selected to determine the spatial pattern of morphologic changes of the landslide. The results show that: (a) A large amount of collapses occurred on dip slopes in the period from 2011 to 2014 and on surrounding debris deposits during the rainy seasons; (b) The average recession rate of the Chunqui Cliff decreased from 24.4 m/yr to 19.6 m/yr compared with the result between 1999 and 2009; (c) The Tsaoling Landslide has lost 6.90 × 106 m3 of soil from November of 2011 to April of 2014, which shows a positive correlation of 0.853 with rainfall; (d) The Chingshui River is undergoing a gradual bed erosion with a volumes of 1.84 × 106 m3.


Introduction
The Tsaoling Landslide, located in central Taiwan, has slid several times at different scales since 1862 [1].One major landslide trigged by the Chi-Chi Earthquake in 1999 has led to a soil loss of about 0.1239 km 3 and notable casualties [2].Moreover, the induced landslide lakes also pose threat to life and property nearby [3].Related researches have been extensively conducted, such as the cause of landslides [4][5][6], volumetric changes [2,7,8], the evolution of topography [1] and the mechanism behind landsides [9,10].However, a series of small-scale landslides occurred after 1999 with significant soil transport was not well studied, which can aslo lead to secondary disasters such as rockslide, dammed lakes and river incision [3,11].Monitoring the spatial pattern and rate of geomorphological changes in this area is indeed vital for protecting residents' life and the ecological environment.
Digital Elevation Model (DEM) or Digital Surface Model (DSM) is a useful tool for the monitoring of large-scale landslides by comparing datasets obtained before and after the events [12][13][14].During the past decades, DEM has been developed from low resolution aerial photo or optical images to high resolution images obtained through airborne LiDAR or Unmanned Aerial Vehicle (UAV), providing more possibilities for monitoring land changes [15,16].However, these measurements have their own limitations, e.g., high cost, sensitive to working conditions or small coverage.On the other hand, InSAR technique, which has been successfully applied to investigate natural hazards such as landslide and earthquakes [17][18][19][20][21][22][23][24], can also be used to generate high quality DEM.Thanks to the notable evolution of sensors and processing techniques, the spatial coverage, resolution, and accuracy of InSAR derived DEM have been steadily improved.Nevertheless DEMs generated by the repeat orbit interferometry still have some constraints, such as spatial-temporal decorrelation, atmospheric artifacts and distortions from other signals (e.g., ground deformation) [25,26].To tackle these constraints, the TanDEM-X (TDX) satellite was launched by DLR and European Aeronautic Defense and Space Company (EADS) Astrium in 2010 to form a constellation with TerraSAR-X (TSX) launched in 2007.As the constellation can acquire an interferometric image pair within several seconds, it is very suitable to generate high quality global DEM [27].A considerable number of studies have been conducted based on TSX/TDX data.For example, 3D topographic reconstruction, fusion of ascending and descending images, phase unwrapping and filtering [28][29][30][31][32], large-scale deformation monitoring e.g., volcano and glacier [33][34][35], relative and absolute accuracy assessment of TSX/TDX DEM [36][37][38][39][40][41], terrain classification and forest height and biomass inversion [42,43].The accumulated TSX/TDX image pairs acquired over the years can be used to generate a set of DEMs over the same area, which is particularly helpful for investigating the temporal evolution of natural hazards, e.g., landslides.To our best knowledge, however, very few researches have been conducted on applications of TSX/TDX DEM time-series in monitoring large-scale landslides.
Focusing on the generation of DEM time series based on TSX/TDX pairs from November of 2011 to April of 2014, this study aims to investigate recent landslides in Tsaoling and to generate the geomorphological changes of landslide.In Section 2, the backgrounds of the Tsaoling landslide is briefly summarized.Section 3 gives an introduction to TSX /TDX datasets and methods as well as an accuracy assessment of the generated DEMs.Results of time-series DEMs in Tsaoling, including volume changes relating to the rainfall and time-series changes of cliffs and Chingshui River, are analyzed in Section 4. Discussions and conclusions are drawn in Sections 5 and 6.

Geological Setting
The study area is located at the joint of Yunlin County and Chiayi County, a foothill region of central Taiwan as outlined by black rectangle in Figure 1.The elevation of this area ranges from 500 m to 1300 m above the sea level (see Figure 1b) and the slopes here commonly dip 30 degrees.The Chingshui River, the only river of this area, flows across the base of slopes and then merges into the Zhuoshui River after it turns north.Taiwan's subtropical monsoon climate brings frequent typhoons (averagely 4 times per year) to this area during the rainy season (May-October) [2].These typhoons are often accompanied by plentiful rainfall (2515 mm/yr), which is another main cause of landslides besides the earthquake [2].The geological formations of Tsaoling, represented by the filled colors on the base map of optical image derived from Google Earth as shown in Figure 2, belonging to three categories: Miocene, Pliocene and Holocene.Miocene has three kinds of terranes including Tawo sandstone, distributing along the Chingshui River with a thickness of 1100 m, Shihliufeng shale and Guandaoshan sandstone.Chingshui shale and Cholan formation are the main terranes of Pliocene and are mainly distributed along the anticline and syncline as shown in Figure 2. Chingshui shale consists of mudstones and shales with intercalated fine-grained sandstone layers.The mudstone has weak cementation and is sensitive to wetting and drying.Cholan formation, which conformably rests on the Chingshui Shale, not only consists of greenish-gray fine-grained sandstones and intercalated shales with high permeability and friability, but also constitutes most of the landslide mass.In addition, landslides before 1979 all occurred in the weak interfaces distributed along Chingshui Shale and Cholan formation [2,8,10].The alluvial deposit, debris deposit and terrace deposit that constitute the latest Holocene with a high compressibility and permeability are the potential localities of collapses and landslides.
Besides the terranes, four faults, one anticline named Tsaoling and one syncline named Chiuhsiungping are also located in the Tsaoling area.The landslides took place in Tsaoling monocline represented by red solid line in Figure 2, dipping towards the south-southwest (SSW).Although the terranes dips towards the SSW on both sides of Chingshui River, the landslides only occurred on the north bank of Chingshui River because the mountain flank dips steeper than the terranes [9].The geological formations of Tsaoling, represented by the filled colors on the base map of optical image derived from Google Earth as shown in Figure 2, belonging to three categories: Miocene, Pliocene and Holocene.Miocene has three kinds of terranes including Tawo sandstone, distributing along the Chingshui River with a thickness of 1100 m, Shihliufeng shale and Guandaoshan sandstone.Chingshui shale and Cholan formation are the main terranes of Pliocene and are mainly distributed along the anticline and syncline as shown in Figure 2. Chingshui shale consists of mudstones and shales with intercalated fine-grained sandstone layers.The mudstone has weak cementation and is sensitive to wetting and drying.Cholan formation, which conformably rests on the Chingshui Shale, not only consists of greenish-gray fine-grained sandstones and intercalated shales with high permeability and friability, but also constitutes most of the landslide mass.In addition, landslides before 1979 all occurred in the weak interfaces distributed along Chingshui Shale and Cholan formation [2,8,10].The alluvial deposit, debris deposit and terrace deposit that constitute the latest Holocene with a high compressibility and permeability are the potential localities of collapses and landslides.
Besides the terranes, four faults, one anticline named Tsaoling and one syncline named Chiuhsiungping are also located in the Tsaoling area.The landslides took place in Tsaoling monocline represented by red solid line in Figure 2, dipping towards the south-southwest (SSW).Although the terranes dips towards the SSW on both sides of Chingshui River, the landslides only occurred on the north bank of Chingshui River because the mountain flank dips steeper than the terranes [9].

TerraSAR-X/TanDEM-X
TerraSAR-X and TanDEM-X satellites having orbits in a closely controlled helix formation act as a large single-pass SAR interferometer to overcome limitations imposed by repeat-pass interferometry (e.g., temporal decorrelation and atmosphere affection) and allow a flexible configuration of the interferometric geometry [27].Three operation modes configured for the systematic data acquisition are the monostatic, alternating monostatic and bistatic modes.The third mode is adopted in this study, as it is commonly used in monitoring land masses, whereas the other two are not designed for DEM generation.It is reported that the absolute vertical height specification accuracy of TSX/TDX DEM is 10 m (90% linear error), and the relative accuracy is ~2 m for areas with predominant terrain slopes that are smaller than 20° [27,37,40,41].

Methodology
Generally speaking, fusion of DEM with different platforms, such as ascending and descending orbits, is a common technique to reconstruct DEM information especially in stable areas with InSAR geometric distortions [28,29].However, the regional terrain in Tsaoling landslide area is too steep (up to 80°) to obtain a complete coverage even using the data from both orbits due to geometry distortions (e.g., layover and shadow).More importantly, the long time span is a vital hamper for reasonable DEM fusion.For instance, the shortest time interval between ascending (28 September 2013) and descending (25 July 2013) orbits is 65 days in this work, spanning a rainy season of Taiwan [2].The heavy rainfall can lead to a huge collapse with a deep of 60 m in a short period [11].Therefore, we have adopted here the traditional differential interferometry of each orbit to analyze the changing features of landslides, where considerable efforts were made with a tradeoff between DEM quality and its spatial resolution.

Methodology and Datasets
3.1.TerraSAR-X/TanDEM-X TerraSAR-X and TanDEM-X satellites having orbits in a closely controlled helix formation act as a large single-pass SAR interferometer to overcome limitations imposed by repeat-pass interferometry (e.g., temporal decorrelation and atmosphere affection) and allow a flexible configuration of the interferometric geometry [27].Three operation modes configured for the systematic data acquisition are the monostatic, alternating monostatic and bistatic modes.The third mode is adopted in this study, as it is commonly used in monitoring land masses, whereas the other two are not designed for DEM generation.It is reported that the absolute vertical height specification accuracy of TSX/TDX DEM is 10 m (90% linear error), and the relative accuracy is ~2 m for areas with predominant terrain slopes that are smaller than 20 • [27,37,40,41].

Methodology
Generally speaking, fusion of DEM with different platforms, such as ascending and descending orbits, is a common technique to reconstruct DEM information especially in stable areas with InSAR geometric distortions [28,29].However, the regional terrain in Tsaoling landslide area is too steep (up to 80 • ) to obtain a complete coverage even using the data from both orbits due to geometry distortions (e.g., layover and shadow).More importantly, the long time span is a vital hamper for reasonable DEM fusion.For instance, the shortest time interval between ascending (28 September 2013) and descending (25 July 2013) orbits is 65 days in this work, spanning a rainy season of Taiwan [2].The heavy rainfall can lead to a huge collapse with a deep of 60 m in a short period [11].Therefore, we have adopted here the traditional differential interferometry of each orbit to analyze the changing features of landslides, where considerable efforts were made with a tradeoff between DEM quality and its spatial resolution.
The flowchart for deriving DEM time series from TSX/TDX image pairs is shown in Figure 3. First, TSX/TDX images with approximate geometrical conditions (i.e., incident angle and pixel resolution) in an independent orbit are coregistered to a reference pair, which is the first acquisition pair in this work.Then, an external coarse-resolution DEM, for example, 1arc s SRTM DEM with a resolution of 30 m, is adopted to remove the topographic phase.As we known, notable resolution difference between external DEM and single-look coregistered slant range single look complex (CoSSC) images (~3 m) will lead to imprecise offset estimation between external DEM and SAR images and localized differential phase error.Hence, a multilooking operation is firstly adopted to improve the coregistration accuracy and reduce phase noise.A non-local filter [44] and a robust coherence calculation [45] are adopted to enhance the success rate of phase unwrapping.Next, the phase unwrapping is directly executed by using minimum cost flow (MCF) method [46].The unwrapped phase is then transferred into heights and a polynomial trend caused by the systematic error is removed from the height map.After this operation, the estimated topographic artefacts are added back to the external DEM to get a refined DEM with higher resolution.Considering that multilooking operation will blur topographic details, an iterative processing that takes the refined DEM as a new external DEM with gradually smaller multilooking factors is suggested to be adopted until the desired spatial resolution of DEM is reached.For a given multilooking factor, the quality of the refined DEM can be assessed with the topographic phase residuals in differential interferograms.Compared to the commonly used multi-looking Differential Interferometric Synthetic Aperture Radar (DInSAR), the iterative operation used here can provide a more accurate DEM with higher resolution.It should be noted that the elevation in these obtained DEMs is with respect to a reference point selected in the phase unwrapping operation, which can be transformed to absolute value with GCPs [37].Finally, a time-series of height differences can be obtained with respect to the selected reference pair.
Remote Sens. 2017, 9, 353 5 of 18 The flowchart for deriving DEM time series from TSX/TDX image pairs is shown in Figure 3. First, TSX/TDX images with approximate geometrical conditions (i.e., incident angle and pixel resolution) in an independent orbit are coregistered to a reference pair, which is the first acquisition pair in this work.Then, an external coarse-resolution DEM, for example, 1arc s SRTM DEM with a resolution of 30 m, is adopted to remove the topographic phase.As we known, notable resolution difference between external DEM and single-look coregistered slant range single look complex (CoSSC) images (~3 m) will lead to imprecise offset estimation between external DEM and SAR images and localized differential phase error.Hence, a multilooking operation is firstly adopted to improve the coregistration accuracy and reduce phase noise.A non-local filter [44] and a robust coherence calculation [45] are adopted to enhance the success rate of phase unwrapping.Next, the phase unwrapping is directly executed by using minimum cost flow (MCF) method [46].The unwrapped phase is then transferred into heights and a polynomial trend caused by the systematic error is removed from the height map.After this operation, the estimated topographic artefacts are added back to the external DEM to get a refined DEM with higher resolution.Considering that multilooking operation will blur topographic details, an iterative processing that takes the refined DEM as a new external DEM with gradually smaller multilooking factors is suggested to be adopted until the desired spatial resolution of DEM is reached.For a given multilooking factor, the quality of the refined DEM can be assessed with the topographic phase residuals in differential interferograms.Compared to the commonly used multi-looking Differential Interferometric Synthetic Aperture Radar (DInSAR), the iterative operation used here can provide a more accurate DEM with higher resolution.It should be noted that the elevation in these obtained DEMs is with respect to a reference point selected in the phase unwrapping operation, which can be transformed to absolute value with GCPs [37].Finally, a time-series of height differences can be obtained with respect to the selected reference pair.

Datasets and Processing
We have collected nine pairs of TSX/TDX data from two orbits, including four pairs from the ascending orbit and five pairs from the descending orbit, as shown in Table 1 and Figure 1.As an example, Figure 4 shows a set of interferograms of descending orbit with different multilooking factors.The areas with low coherence (<0.18) and serious geometric distortion were masked out aiming to improve the quality of phase unwrapping.Notable topographic phase residuals can be seen in Figure 4c, which is because the topography was altered by landslides occurred in the period

Datasets and Processing
We have collected nine pairs of TSX/TDX data from two orbits, including four pairs from the ascending orbit and five pairs from the descending orbit, as shown in Table 1 and Figure 1.As an example, Figure 4 shows a set of interferograms of descending orbit with different multilooking factors.The areas with low coherence (<0.18) and serious geometric distortion were masked out aiming to improve the quality of phase unwrapping.Notable topographic phase residuals can be seen in Figure 4c, which is because the topography was altered by landslides occurred in the period of 2000 and 2011 [3,7,11].In the first iterative full resolution differential interferogram (i.e., Figure 4d), the topographic residual is still visible, which is mainly due to resolution discrepancy between updated external DEM (6 m) and the interferogram (3 m).By further updating the external DEM according to the phase residuals in the first iteration, the topographic residuals in the second iterative interferogram (i.e., Figure 4e) become ignorable, indicating that after the first iteration, the topography over interest area has already been fully reconstructed.Similarly, the time-series DEMs of ascending orbit is also obtained after 1st iteration.In addition to SAR images, we have also collected rainfall datasets between November 2011 and May 2014 from the Tongtou Station, which is about 6 km away from Tsaoling Landslide.The histogram is shown in Figure 5.The blue peaks represent typhoons that had visited Taiwan (8 times in 2012 and 6 times in 2013).All of the typhoons occurred during the period from June to October and were accompanied by abundant rainfall. of 2000 and 2011 [3,7,11].In the first iterative full resolution differential interferogram (i.e., Figure 4d), the topographic residual is still visible, which is mainly due to resolution discrepancy between updated external DEM (6 m) and the interferogram (3 m).By further updating the external DEM according to the phase residuals in the first iteration, the topographic residuals in the second iterative interferogram (i.e., Figure 4e) become ignorable, indicating that after the first iteration, the topography over interest area has already been fully reconstructed.Similarly, the time-series DEMs of ascending orbit is also obtained after 1st iteration.In addition to SAR images, we have also collected rainfall datasets between November 2011 and May 2014 from the Tongtou Station, which is about 6 km away from Tsaoling Landslide.The histogram is shown in Figure 5.The blue peaks represent typhoons that had visited Taiwan (8 times in 2012 and 6 times in 2013).All of the typhoons occurred during the period from June to October and were accompanied by abundant rainfall.

Relative Height Accuracy
Absolute height accuracy is a common criterion to assess InSAR DEMs based on the accurate traditional outcomes, such as GPS and LiDAR [40,41].However, in some areas, ground control points is lacking or difficult to collect.Hence, the relative accuracy of DEM achieved from repeated observations or high quality points in the study area may be an important quality indicator [36][37][38].The main procedure of relative accuracy calculation in this study is as follows.We firstly chose 27 uniformly distributed stable points together with their height outside the active landslide area, based on the previous work [11], slope maps as well as coherence maps of Tsaoling Landslide.A total of 16 (out of 27) points represented by the black filled rectangles in Figure 2 are used to mitigate the systematic errors in InSAR derived DEMs with polynomials described in Gruber et al. [37].The remained 11 marked with black filled triangles in Figure 2 are used for validation by calculating the RMSE of the height difference between calibrated values.As images obtained on 11 November 2013 and 16 January 2014 are not sufficient to cover the entire interest area, the accuracy of these two pairs are calculated based on 5 stable points in overlapped area.As shown in Table 2, the relative height accuracy of the descending orbit are 0.98 m, 1.27 m, 1.12 m, 1.00 m and 1.90 m with an average value of 1.30 m, and the relative values of the ascending orbit are 1.14 m, 1.44 m, 1.84 m and 2.25 m with an average value of 1.12 m.

Characteristics of DEM Time Series
Seven InSAR derived DEMs of the Tsaoling Landslide are showed in Figure 6. Figure 6a shows the difference in DEMs collected in December of 2012 and September of 2013.There are some voids due to the use of single pass side-looking data, while it is evident that erosions are mainly distributed in debris deposit (Figure 2) and are concentrated on two banks of the Chingshui River, where two

Relative Height Accuracy
Absolute height accuracy is a common criterion to assess InSAR DEMs based on the accurate traditional outcomes, such as GPS and LiDAR [40,41].However, in some areas, ground control points is lacking or difficult to collect.Hence, the relative accuracy of DEM achieved from repeated observations or high quality points in the study area may be an important quality indicator [36][37][38].The main procedure of relative accuracy calculation in this study is as follows.We firstly chose 27 uniformly distributed stable points together with their height outside the active landslide area, based on the previous work [11], slope maps as well as coherence maps of Tsaoling Landslide.A total of 16 (out of 27) points represented by the black filled rectangles in Figure 2 are used to mitigate the systematic errors in InSAR derived DEMs with polynomials described in Gruber et al. [37].The remained 11 marked with black filled triangles in Figure 2   Seven InSAR derived DEMs of the Tsaoling Landslide are showed in Figure 6. Figure 6a shows the difference in DEMs collected in December of 2012 and September of 2013.There are some voids due to the use of single pass side-looking data, while it is evident that erosions are mainly distributed in debris deposit (Figure 2) and are concentrated on two banks of the Chingshui River, where two slope failures P and Q are observable.The deepest collapse depth of P is about 50 m.Besides the erosion along the Chingshui River, two rockslides located on the Chunqui cliff (N) and cliff J are observed on the dip slope.Figure 6b,c    The landslide area is bounded by closed black solid lines in Figure 6.Two profiles marked in black solid lines M and N, and four active slope failures marked in blue signs P, Q, I, and J in Figure 2 were chosen to conduct the analysis.
Figure 6d-g are results of the descending orbit, indicating that voids located on the active area will lead to an unreliable DEM fusion with two or more platforms.Figure 6d,e have shown that the erosion is mainly located on the Chunqui Cliff, resting on the Chingshui shale, whose magnitude is generally smaller than 15 m and can reach 60 m at maximum.The depth of deposit for Chunqui cliff is smaller than 8 m with certain areas deeper than 30 m.In addition to the Chunqui Cliff (N), there is a small-scale collapse on cliff J with an erosion of 20 m at maximum and two obvious erosion areas outside the interest area that can be recognized in optical images.For example, one small-scale collapse represented by L1 with slip depth less than 26 m (Figure 6g,h) and another artificial topographic changing is observed and marked by L2 (Figure 6g,i).Figure 6f,g show the erosion on five selected cliffs, N, P, Q, I and J.The maximum depth of erosion on the Chunqui Cliff (N) is 74 m and the collapse trend resembles the Cliff's strike direction.Cliffs P, Q and I, distributed along the Chingshui River, have obvious erosion, especially cliff P.However, the topographic changes of cliff I was incomplete due to the geometric distortion of InSAR technique.Besides the five cliffs, obvious erosion was also observed in the channel located opposite cliff I.The depth of erosion in the channel is less than 15 m except for partially areas whose depth is as large as 32 m.In addition to the erosions, evident deposits are observed below the Chunqui Cliff (N) with a maximum height of 35 m.

Analysis of Adjacent DEMs
This study also obtained the height difference between adjacent pairs, as showed in Figure 7, to get detailed features of collapses.Three adjacent intervals of ascending orbit are shown in Figure 7a-c  Figure 6d-g are results of the descending orbit, indicating that voids located on the active area will lead to an unreliable DEM fusion with two or more platforms.Figure 6d,e have shown that the erosion is mainly located on the Chunqui Cliff, resting on the Chingshui shale, whose magnitude is generally smaller than 15 m and can reach 60 m at maximum.The depth of deposit for Chunqui cliff is smaller than 8 m with certain areas deeper than 30 m.In addition to the Chunqui Cliff (N), there is a small-scale collapse on cliff J with an erosion of 20 m at maximum and two obvious erosion areas outside the interest area that can be recognized in optical images.For example, one small-scale collapse represented by L1 with slip depth less than 26 m (Figure 6g,h) and another artificial topographic changing is observed and marked by L2 (Figure 6g,i).Figure 6f,g show the erosion on five selected cliffs, N, P, Q, I and J.The maximum depth of erosion on the Chunqui Cliff (N) is 74 m and the collapse trend resembles the Cliff's strike direction.Cliffs P, Q and I, distributed along the Chingshui River, have obvious erosion, especially cliff P.However, the topographic changes of cliff I was incomplete due to the geometric distortion of InSAR technique.Besides the five cliffs, obvious erosion was also observed in the channel located opposite cliff I.The depth of erosion in the channel is less than 15 m except for partially areas whose depth is as large as 32 m.In addition to the erosions, evident deposits are observed below the Chunqui Cliff (N) with a maximum height of 35 m.

Analysis of Adjacent DEMs
This study also obtained the height difference between adjacent pairs, as showed in Figure 7, to get detailed features of collapses.Three adjacent intervals of ascending orbit are shown in Figure 7a-c and four intervals of the descending orbit showed in Figure 7d-g.Figure 7a spans the whole rainy season, during which five typhoons occurred with abundant rainfall (the average daily precipitation is about 284 mm; See Figure 4).Two major collapses (cliff P with maximum depth of up to 50 m and cliff Q) are observed.Small-scale collapses are observed on cliff J and the Chunqui Cliff (N), whose maximum collapse depth has reached 40 m.In addition, river incision is observed in the upstream with a maximum dissection of 18 m.However, there is no obvious collapse on the slope or along the river because the data of Figure 7b,c were collected in dry season.
Data of Figure 7d,e were collected in dry season from November of 2011 to March of 2012, with a total rainfall of 35.5 mm and 148.5 mm respectively.Figure 7d,e indicate that no evident erosions and deposition are observed except a new branch of Chunqui Cliff (N), represented by red dotted lines in Figure 7d,e, and some small-scale erosions on cliff J.
The deformation of the interested area mainly occurred during March of 2012 to July of 2013 (Figure 7f), when ten typhoons came with heavy rains and the maximum daily rainfall was about 565 mm.Collapses are mainly located on cliff N, P, Q, I and J, where the maximum depth of collapse is 60 m located on cliff N. Lots of small-scale erosions situated between the Chunqui Cliff (N) and the Chingshui River are also observed whose erosion depth is less than 12 m.Four erosion areas along the river are observed and are represented by black dotted lines in Figure 7f.Three of them are located in the upstream and one in the downstream.Besides collapses, obvious deposition distributed along the cliffline of the Chunqui Cliff with a maximum depth of 30 m are also observed.
During July of 2013 to April of 2014, four recorded typhoons brought a total rainfall of 1961 mm.Erosion and deposition are mainly occurred on cliff P, Q and J and their underneath in Figure 7g, besides some small-scale changes in cliff I, indicating the instability of these cliffs.Meanwhile, smaller erosions in the Chingshui River and along its banks are observed.More details will be depicted in Section 4.1.3.

Changes of the Chingshui Channel
In addition to the analysis of dip slope, the Chingshui River flowing across the interested area, represented by cyan filled area in Figure 2, is extracted to analyze detailed topographic changes in Figure 8. Height difference among DEMs generated from adjacent pairs on the ascending orbit are shown in Figure 8a-c and four on the descending orbit are shown in Figure 8d-g.The incision of the river mainly occurred during March of 2012 to July of 2013 (Figure 8f).Three sections incised over 20 m are marked by red solid lines in Figure 8f, two of which are located at corners of the river.At those corners, the river incision rate and strength are huge and resulted in a collapse of riverbank.The rest one is located at the middle part of these two corners.Furthermore, river incision created during July of 2013 to April of 2014 (Figure 8g) is smaller than that created during March of 2012 to July of 2013 (Figure 8f) and is mainly concentrated in the middle of the river instead of the banks.
Erosions of the Chingshui River on the ascending orbit are shown in Figure 8a-c.An evident erosion occurred during December of 2012 to September of 2013 (Figure 8a), and the strength of incision at corners of the Chingshui River is far less, indicating that the incision in the downstream happened during March and December of 2012.The remaining four periods (Figure 8b-e) that belong to the dry season have rather stable conditions apart from some small-scale local erosions and river incisions.

Analysis of Two Selected Profiles
To figure out the changing features of cliffs of the Tsaoling Landslide, two profiles (M and N) across the Chunqui Cliff (N) and cliff Q in Figure 7    Two bowls, f1 and f2, on cliff P and cliff N are observed in Figure 9a,c.The variable DEM difference of f1 indicates that cliff P is still active, and the aquamarine filled space refers to continual erosion and the yellow filled space refers to a second collapse under the previous one.Similarly, f2 is experiencing erosion and deposit.However, a new bowl f3 in Figure 9c cropped up during March and December of 2012 and the erosion is continuing.On the other hand, the elevation difference of profile N increases with time as shown in Figure 9b,d, confirming that the Chunqui Cliff (N) is undergoing erosions especially during March and July of 2013.

Change Rate of Cliffs
To quantify the change rate of cliffs in the Tsaoling Landslide, the recession area and rate of the cliffs are analyzed according to five TSX/TDX pairs on the descending orbit due to the full coverage over the interest area through the time.As described in Section 3.2, the first pair acquired at November 2011 is selected as the reference pair.It is noting that the erosion boundaries of these cliffs were created by manual segmentation with a threshold of 8 m in height difference maps.We do not adopted the visual interpretation here, similar to optical images [7,11], to obtain the boundaries in the intensity maps, due to the confusion caused by geometric distortion in side-looking geometry.Figure 10 shows the history of cliffs' boundaries with an optical image as the background.On the other hand, cliff P and Q do not have four observations because they appeared after March 2012.It is evident that erosions of these five cliffs mainly occurred during November of 2011 to July of 2013, which are represented by black solid lines in Figure 10.Cliff N and I also have small-scale erosions during July of 2013 to April of 2014, whereas cliff P experienced a large-scale erosion together with cliff Q on the northeast in the same period as marked green solid lines in Figure 10.
Two bowls, f1 and f2, on cliff P and cliff N are observed in Figure 9a,c.The variable DEM difference of f1 indicates that cliff P is still active, and the aquamarine filled space refers to continual erosion and the yellow filled space refers to a second collapse under the previous one.Similarly, f2 is experiencing erosion and deposit.However, a new bowl f3 in Figure 9c cropped up during March and December of 2012 and the erosion is continuing.On the other hand, the elevation difference of profile N increases with time as shown in Figure 9b,d, confirming that the Chunqui Cliff (N) is undergoing erosions especially during March and July of 2013.

Change Rate of Cliffs
To quantify the change rate of cliffs in the Tsaoling Landslide, the recession area and rate of the cliffs are analyzed according to five TSX/TDX pairs on the descending orbit due to the full coverage over the interest area through the time.As described in Section 3.2, the first pair acquired at November 2011 is selected as the reference pair.It is noting that the erosion boundaries of these cliffs were created by manual segmentation with a threshold of 8 m in height difference maps.We do not adopted the visual interpretation here, similar to optical images [7,11], to obtain the boundaries in the intensity maps, due to the confusion caused by geometric distortion in side-looking geometry.Figure 10 shows the history of cliffs' boundaries with an optical image as the background.On the other hand, cliff P and Q do not have four observations because they appeared after March 2012.It is evident that erosions of these five cliffs mainly occurred during November of 2011 to July of 2013, which are represented by black solid lines in Figure 10.Cliff N and I also have small-scale erosions during July of 2013 to April of 2014, whereas cliff P experienced a large-scale erosion together with cliff Q on the northeast in the same period as marked green solid lines in Figure 10.The recession area of each cliff is calculated with the reference period between November of 2011 and December of 2011 (represented by red solid lines in Figure 10) instead of the boundaries in reference date, i.e., December of 2011, due to the unrecognizable boundaries in intensity map aforementioned.The recession areas of cliff P and Q are calculated based on pairs collected between November of 2011 and July of 2013 (represented by black solid lines in Figure 10) because of the stability before March of 2012.The history of recession areas are calculated according to the pixel size as shown in Table 3.The erosion of cliff N and I occurred during March of 2012 to July of 2013 (Figure 7f) are sized 31.22 × 10 3 m 2 and 14.16 × 10 3 m 2 .The erosion of Cliff P and Q occurred after March 2012 and the changing areas are 31.00× 10 3 m 2 and 7.06 × 10 3 m 2 .However, cliff J has undergone two stages: the first one is the increasing erosion area from 2.49 × 10 3 m 2 to 4.63 ×10 3 m 2 at the top of cliff J, and the other one is a secondary erosion rested on the previous one with an area of 4.06 ×10 3 m 2 .Besides the size, changing rates of these areas are also calculated as shown in Table 3.Two of five cliffs retreated with a decreasing rate, for example, the recession rate of cliff N decreased from 38.14 × 10 3 m 2 /yr to 15.43 ×10 3 m 2 /yr, and cliff I from 23.72 ×10 3 m 2 /y to 11.56 ×10 3 m 2 /yr.On the other hand, cliff P retreats at a rate of 44.72 ×10 3 m 2 /yr, whereas cliff Q has a relatively lower rate about 10.18 ×10 3 m 2 /yr, indicating its stability.In addition to the recession area and rate, the maximum retreat rate of each cliff which is defined along the maximum recession direction of the cliff [7], represented by S1 to S5 in Figure 10, is estimated.The maximum retreat rate can be regarded as an assessment indicator of the whole cliff's recession [7,11].Table 4 shows the retreat distance and rate of these cliffs.The retreat rates of cliff N and I reached their peaks of 26.6 m/yr and 19 m/yr, respectively, during March of 2012 to July of 2013.Besides, the retreat rates of cliff P (111.43 m/yr) and Q (41.47 m/yr) are relatively high during period July of 2013 to April of 2014 (Figure 7g).Cliff J retreated with a distance of 107.3 m during period December of 2011 to July of 2013, and then a secondary collapse occurred with a distance of 85.01 m.The average rates of these cliffs (N, P, Q, and I) are 19.73 m/yr, 111.43 m/yr, 41.47 m/yr and 12.76 m/yr.However, the rate of cliff J is not calculated, due to its specificity aforementioned.Compared to the results of 2009 [7,11], the retreat rate of the Chunqui Cliff (N) has decreased from 24.4 m/yr to 19.73 m/yr.In addition, the maximum retreat distance of cliff J moved from the northwest to the southeast during November of 2011 to April of 2014 see the red arrow in Figure 10b.

Calculation of the Changing Volumes
Traditional DEMs derived from different platforms, such as LiDAR and optical images, need to take the coregistrated or oversampled error into consideration because of the inconsistent resolutions.In contrast, InSAR derived DEMs are coregistered with sub pixel accuracy, resulting in a high precision of relative volumetric assessment.Volumetric changes of soil in the Tsaoling Landslide on the two orbits (ascending and descending) are calculated by the multiplication between pixel area, i.e., ground range resolution × azimuth resolution, and its height difference and shown in Table 5.The values in this table are the sum total volumes of pixels in TSX/TDX pairs, except for the voids caused by geometric distortion.The results show that the maximum volume changes on the descending orbit occurred during March of 2012 and July of 2013 (Figure 7f) with an erosion of 10.23 ×10 6 m 3 and a deposit of 6.57 ×10 6 m 3 .On the other hand, the minimum change happened during November of 2011 to December of 2011 (Figure 7d) with an erosion of 5.58 ×10 6 m 3 and deposit of 4.42 ×10 6 m 3 .The changed volumes on the descending orbit were −1.16 ×10 6 m 3 , 0.40 ×10 6 m 3 , −3.66 ×10 6 m 3 and −2.48 ×10 6 m 3 as shown in Figure 7d-g with a total of −6.90 ×10 6 m 3 .The −2.20 ×10 6 m 3 on the ascending orbit indicates that the soil of landslide is in a state of loss.Moreover, a regression analysis is conducted to obtain the relationship between rainfall and volumes.The correlation coefficient among the erosion volumes, deposit volumes and changing volumes are 0.868, 0.349 and 0.853 respectively, indicating a positive correlation between soil erosion and rainfall.The low coefficient between rainfall and deposit volumes may be caused by the transport of Chingshui river.

Continual Incision of the Chingshui River
The stability of the Chingshui River is an important safety indicator of water and soil, as well as the life and property around it [3].Therefore, the volume history of the Chingshui River during the acquisition time are analyzed independently, which can also be defined as the river's transport capacity due to the stability of the river's left bank as described in Section 2 [9].The results are shown in the 5th column of Table 6, which should also equal to the volume difference between the whole interest area and Chingshui River considering that the river is the only way to transport the soil in this area.It is clear that the maximum transport capability is observed during July of 2013 to April of 2014 (−3.26 ×10 6 m 3 , Figure 8g) and the minimum one during November of 2011 to December of 2011 (1.28 ×10 6 m 3 , Figure 8d).Moreover, a small-scale deposit happened during December of 2011 to March of 2012 in the downstream in Figure 8e is 1.67 ×10 6 m 3 .In addition, the maximum incision occurred during March of 2012 to July of 2013, is 1.47 ×10 6 m 3 , as shown in Figure 8f.The 4th column of Table 6 indicate that the Chingshui River incised with a total volume of 1.84 ×10 6 m 3 between November of 2011 and April of 2014, where the calculation is based on the data acquired along the descending orbit.

Discussion
In general, lots of studies have been done to obtain the dynamic activities i.e., volume changing, of Tsaoling landslide in the past decade based on DEMs derived from LiDAR or optical images [2,3,7,8,11].Unlike these methods, this study used an iterated differential interferometry SAR method described in Section 3 to generate TSX/TDX DEM time-series with high spatial resolution (~3 m) and relative precision (~1.2 m) over the period from November of 2011 to April of 2014.Furthermore, compared with the existed LiDAR DEMs (time interval larger than 1 year) [2,7,11], the TSX/TDX DEM time-series can generate more details about geomorphological changes of Tsaoling landslide because of its much smaller time interval (the minimum interval equals to 22 days).This study provides case example to evaluate the impacts of landslides on water source planning and hazard assessment as well as the importance of geological factors in landslide occurrence.
In our study, a tradeoff between resolution and geometric distortion should be taken into consideration.As we known, the effect of geometric distortion in differential interferograms is increasing with the DEM resolution through an iterated operation [28], resulting in voids for volume calculation of erosion and deopition in landslide areas.However, detail information will be lost if we do not adopted the iteration.Fortunately, geometric distortion of Tsaoling landslide is not significant due to the smooth of the slide surface except for few steep cliffs.Furthermore, SAR images from different orbits (ascending and descending) are collected to fill the gaps caused by geometric distortion.Therefore, the iterated operation is adopted to obtain the final DEMs.The volume history for each orbit is calculated without the values in voids areas, small areas whose erosion or deposition are not obvious in Figure 7.
Like the landslides occurred in other areas [18], we confirmed that slope failures during the period from November of 2011 to April of 2014 in Tsaoling landslide area are mainly triggered by events, such as typhoons, accompanied with abundant rainfall.The results suggest a positive correlation (coefficient equals to 0.868) between erosion volume and rainfall, and the slope failures are mainly occurred during the rainy season from March of 2012 to July of 2013 with a maximum collapse depth of 70 m.
Regarding the recession of Chunqui cliff after 1999 Chi-Chi earthquake, our conclusions are somewhat different.The comparison between the average recession rate calculated in this study (19.73 m/yr) and the corresponding rate (24.4 m/yr) as determined by Kuo [7] during the period from 1999 to 2009, indicates that a decreasing recession of Chunqui cliff.This result is in accordance with the prediction of exponential regression [7].Besides Chunqui cliff, the recession rates of other four selected (i.e., cliff P, Q, I and J) cliffs are 111.43, 41.47, 12.76 and 65.93 m/yr, indicating the instabilities in the channel of Chingshui River and cliff J [3].We also calculated the erosion and deposition volumes of the study area, 28.43 ×10 6 m 3 for erosion and 21.53 ×10 6 m 3 for deposition, the volume difference equals to 6.90 ×10 6 m 3 , translating to annual erosion rate equal to 2.90 ×10 6 m 3 /yr.Compared with the results between 2000 and 2003, i.e., 10 ×10 6 m 3 /yr [2], the erosion rate is also decreasing obviously.This is probably a consequence of a decreasing deposition by transportation of Chingshui River within 15 years.Moreover, the volume difference of Chingshui River is also calculated, i.e., 1.84 ×10 6 m 3 , concentrating on the upstream with a maximum erosion depth of 20 m.The erosion depth is larger than the predicted results obtained during period from 2001 to 2005 (i.e., 17 m) [3].Besides the erosion, small-scale deposition is also observed during the dry season without great events, see Table 6.The volumes difference indicates that the Chingshui River section selected in our study area is undergoing an alternating erosion-deposition pattern, as predicted in the study of Hsu [3].

Conclusions
This paper used an iterated DInSAR techniques to generate time-series DEMs of Tsaoling landslide with two different platforms TSX/TDX SAR images over the period from November of 2011 to April of 2014.Geomorphological changes and the current state of the landslide area are analyzed.The results have shown that erosions and deposits mainly occurred on the unstable geological stratals, such as Debris deposit, Cholan formation and Chingshui shale, during the rainy season when typhoons stroke the area with abundant rainfall.Five cliffs and the Chingshui River were selected to analyze the changing feature of terrain in this area.We have found that the largest erosion and collapse, which is about 70 m, occurred during March of 2012 to July of 2013 accompanying with small-scale collapses.The retreat rate of the most representative cliff, the Chunqui Cliff (N), decreased from 24.4 m/yr to 19.6 m/yr compared with the data collected between 1999 and 2009, whereas the retreat distance of cliff J continued to decrease until the secondary collapse happened in the same area.In addition, the landslide mass in the study area continued to losing according to the estimated volume change, i.e., 6.90 ×10 6 m 3 for descending platform and 2.20 ×10 6 m 3 for ascending platform, which is positively correlated with rainfall with a coefficient of 0.853.Moreover, the Chingshui River is undergoing a gradual bed erosion with a total volume of 1.84 ×10 6 m 3 .The encouraging results indicate that TSX/TDX DEM time-series might be a promising choice for tracking the movement of landslides, especially for rapid one.In the next step, time-series DEMs with large coverage will be obtained to analyze the transportation of landslide mass widely.

Figure 1 .
Figure 1.Location and data coverage of the Tsaoling Landslide.Blue rectangles and red rectangles in (a) represent the descending and ascending TSX/TDX pairs respectively.The black rectangle is the location of the Tsaoling Landslide with an enlarged version as shown in (b).Blue rectangles in (b).They are the clipped area of two orbits and the black rectangle represents the final interest area.

Figure 1 .
Figure 1.Location and data coverage of the Tsaoling Landslide.Blue rectangles and red rectangles in (a) represent the descending and ascending TSX/TDX pairs respectively.The black rectangle is the location of the Tsaoling Landslide with an enlarged version as shown in (b).Blue rectangles in (b).They are the clipped area of two orbits and the black rectangle represents the final interest area.

Figure 2 .
Figure 2. Geology map of the investigated area (modified from Chen et al. 2006 and Tang et al. 2009).

Figure 4 .
Figure 4. DEM reconstruction by iterative DInSAR.(a) Original interferogram from image pair acquired on 17 November 2011 where 2 × 2 multilooking operation is performed; (b) Flattened interferogram; (c) Original DInSAR result; (d) 1st full resolution iterative DInSAR result with the updated external DEM; (e) 2nd full resolution iterative DInSAR result where the topographic artefacts become ignorable; (f) full resolution coherence map.The background of (a-f) is the average intensity map of TSX/TDX.The black holes in (c-e) are the masked areas with low coherence and serious geometric distortion.

Figure 4 .
Figure 4. DEM reconstruction by iterative DInSAR.(a) Original interferogram from image pair acquired on 17 November 2011 where 2 × 2 multilooking operation is performed; (b) Flattened interferogram; (c) Original DInSAR result; (d) 1st full resolution iterative DInSAR result with the updated external DEM; (e) 2nd full resolution iterative DInSAR result where the topographic artefacts become ignorable; (f) full resolution coherence map.The background of (a-f) is the average intensity map of TSX/TDX.The black holes in (c-e) are the masked areas with low coherence and serious geometric distortion.

Figure 5 .
Figure 5. Histogram of daily rainfall (observed at the Tongtou station close to the Tsaoling Landslide).

Figure 5 .
Figure 5. Histogram of daily rainfall (observed at the Tongtou station close to the Tsaoling Landslide).
are used for validation by calculating the RMSE of the height difference between calibrated values.As images obtained on 11 November 2013 and 16 January 2014 are not sufficient to cover the entire interest area, the accuracy of these two pairs are calculated based on 5 stable points in overlapped area.As shown in Table2, the relative height accuracy of the descending orbit are 0.98 m, 1.27 m, 1.12 m, 1.00 m and 1.90 m with an average value of 1.30 m, and the relative values of the ascending orbit are 1.14 m, 1.44 m, 1.84 m and 2.25 m with an average value of 1.12 m.
show a difference in height of DEM collected in November of 2013 and January of 2014.Deposit and erosion were found mainly around the Chingshui River.slope failures P and Q are observable.The deepest collapse depth of P is about 50 m.Besides the erosion along the Chingshui River, two rockslides located on the Chunqui cliff (N) and cliff J are observed on the dip slope.Figure 6b,c show a difference in height of DEM collected in November of 2013 and January of 2014.Deposit and erosion were found mainly around the Chingshui River.

Figure 6 .
Figure 6.DEM time series from the ascending and descending orbits.The accumulated area is represented by yellow and erosion area is represented by purple.(a-c): the elevation difference on the ascending orbit based on a reference date of 16 December 2012;(d-g): the elevation difference on the descending orbit based on a reference date of 17 November 2011.(h,i) are two selected areas with obvious elevation change outside the landslide area as showed with red solid lines in (e-g).The landslide area is bounded by closed black solid lines in Figure6.Two profiles marked in black solid lines M and N, and four active slope failures marked in blue signs P, Q, I, and J in Figure2were chosen to conduct the analysis.

Figure 6 .
Figure 6.DEM time series from the ascending and descending orbits.The accumulated area is represented by yellow and erosion area is represented by purple.(a-c): the elevation difference on the ascending orbit based on a reference date of 16 December 2012; (d-g): the elevation difference on the descending orbit based on a reference date of 17 November 2011.(h,i) are two selected areas with obvious elevation change outside the landslide area as showed with red solid lines in (e-g).The landslide area is bounded by closed black solid lines in Figure6.Two profiles marked in black solid lines M and N, and four active slope failures marked in blue signs P, Q, I, and J in Figure2were chosen to conduct the analysis.
and four intervals of the descending orbit showed in Figure7d-g.Remote Sens. 2017, 9, 353 9 of 18

Figure 7 .
Figure 7. Elevation difference between two adjacent dates for two orbits.(a-c): elevation difference from December of 2012 to September of 2013, September of 2013 to November of 2013 and November of 2013 to January of 2014 respectively for the ascending orbit of investigated area; (d-g): elevation difference from November of 2011 to December of 2011, December of 2011 to March of 2012, March of 2012 to July of 2013 and July of 2013 to April of 2014 respectively for the descending orbit.

Figure 7 .
Figure 7. Elevation difference between two adjacent dates for two orbits.(a-c): elevation difference from December of 2012 to September of 2013, September of 2013 to November of 2013 and November of 2013 to January of 2014 respectively for the ascending orbit of investigated area; (d-g): elevation difference from November of 2011 to December of 2011, December of 2011 to March of 2012, March of 2012 to July of 2013 and July of 2013 to April of 2014 respectively for the descending orbit.

Figure 8 .
Figure 8. Elevation difference of adjacent TSX/TDX pairs of the Chingshui River.(a-c): on the ascending orbit; (d-g): on the descending orbit.The red solid lines represent the obvious incisions of Chingshui River.

4. 2 . 1 .
Analysis of Two Selected ProfilesTo figure out the changing features of cliffs of the Tsaoling Landslide, two profiles (M and N) across the Chunqui Cliff (N) and cliff Q in Figure7are selected.Profile M is similar to the one commonly used by researchers and it passes through the deepest part of cliff Q and cliff N [2,9].On the other hand, profile N distributed along the direction of the Chunqui cliffline, helping to quantitatively analyze the changing history.Figure 9a,c where the reference pairs were collected in December of 2012 and November of 2011 respectively indicate the results of profile M for the ascending and descending orbits.Figure 9b,d represent the results of profile N for the two orbits with the same reference pairs.

Figure 9 .
Figure 9. DEM difference of profile M (a,c) and profile N (b,d) in Figure 6 for ascending and descending orbits respectively.Red, blue, green solid lines in (a,b) represent the results of Figure 6a-c respectively, whereas red, blue, black and green solid lines in (c,d) represent the results of Figure 6d-g.

Figure 8 .
Figure 8. Elevation difference of adjacent TSX/TDX pairs of the Chingshui River.(a-c): on the ascending orbit; (d-g): on the descending orbit.The red solid lines represent the obvious incisions of Chingshui River.

4. 2 . 18 Figure 8 .
Figure 8. Elevation difference of adjacent TSX/TDX pairs of the Chingshui River.(a-c): on the ascending orbit; (d-g): on the descending orbit.The red solid lines represent the obvious incisions of Chingshui River.
are selected.Profile M is similar to the one commonly used by researchers and it passes through the deepest part of cliff Q and cliff N [2,9].On the other hand, profile N distributed along the direction of the Chunqui cliffline, helping to quantitatively analyze the changing history.Figure 9a,c where the reference pairs were collected in December of 2012 and November of 2011 respectively indicate the results of profile M for the ascending and descending orbits.Figure 9b,d represent the results of profile N for the two orbits with the same reference pairs.

Figure 9 .
Figure 9. DEM difference of profile M (a,c) and profile N (b,d) in Figure 6 for ascending and descending orbits respectively.Red, blue, green solid lines in (a,b) represent the results of Figure 6a-c respectively, whereas red, blue, black and green solid lines in (c,d) represent the results of Figure 6d-g.

Figure 9 .
Figure 9. DEM difference of profile M (a,c) and profile N (b,d) in Figure 6 for ascending and descending orbits respectively.Red, blue, green solid lines in (a,b) represent the results of Figure 6a-c respectively, whereas red, blue, black and green solid lines in (c,d) represent the results of Figure 6d-g.

Figure 10 .
Figure 10.Recession of the erosion boundaries of five selected cliffs.Red, blue, black and green lines represent the erosion boundaries of Figure 6d-g respectively.The base map of optical image derived from Google Earth (28 November 2015), whose coverage is represented by two black rectangles in Figure 2.

Figure 10 .
Figure 10.Recession of the erosion boundaries of five selected cliffs.Red, blue, black and green lines represent the erosion boundaries of Figure 6d-g respectively.The base map of optical image derived from Google Earth (28 November 2015), whose coverage is represented by two black rectangles in Figure 2.

Table 2 .
The relative accuracy of TSX/TDX DEMs.

Table 2 .
The relative accuracy of TSX/TDX DEMs.

Table 3 .
Recession area and rate of five selected cliffs.Positive values means the retreat of cliffs whereas negative values means the opposite direction of retreat.

Table 4 .
Retreat distance and rate of each cliffline.
Positive values means the retreat of cliffs whereas negative values means the opposite direction of retreat.

Table 5 .
Volumetric change and rainfall in the investigated area.
Note: The negative values in column 4 are referred to the erosion, whereas the positive value means deposit.

Table 6 .
Volume change of the Chingshui River.