Investigation of Slow-Moving Artificial Slope Failure with Multi-Temporal InSAR by Combining Persistent and Distributed Scatterers: A Case Study in Northern Taiwan

In Taiwan, landslides pose serious threats to local residents and infrastructures each year due to high mountain relief and distinct seasonal precipitation distribution. Interferometric synthetic aperture (InSAR) provides a powerful tool to map terrain motion and characterize the failure mechanism of unstable slopes. However, it is challenging for the conventional InSAR technique to obtain reliable landslide information in mountainous regions because of insufficient coherent measurements and signal confusion caused by vegetation coverage and rugged terrain. In this study, we adopt an optimized multi-temporal InSAR (MTInSAR) approach to analyze the surface displacement of an artificial side slope along Freeway No. 3, where a catastrophic landslide failure occurred on 25 April 2010, in northern Taiwan. To increase the spatial extent of the deformation signal, we integrate information from both persistent scatterers (PSs) and distributed scatterers (DSs). Topographic residual and height-dependent atmospheric delays are corrected by a component-based method and joint model estimation, respectively. The results reveal the existence of slope movement with a rate of about −30 mm/year prior to the landslide failure. Further analysis shows that the temporal behaviors of downslope movement are correlated with local precipitation. The study demonstrates the need to continuously monitor and verify the stability of artificial slopes to prevent and minimize the probability of a similar landslide occurrence in the future.


Introduction
Landslides, as the result of gravitational force and geological instability, are one of the major types of geological disasters around the world [1,2]. To mitigate and reduce landslide hazards, enormous effort has been dedicated to investigating the failure mechanism and assessing the potential vulnerability, such as landslide inventories, forensic assessment, and displacement inventories [3][4][5][6]. However, the current forecasting of landslide location and timing remains a challenge, primarily due to the complex relationship between the abrupt landslide movement and underlying triggers (e.g., precipitation history and seismic events) [7]. Alternatively, landslide deformation mapping over time provides critical information for the risk assessment of the impending slope failure [7,8].
Conventional geodetic techniques (e.g., Global Positioning System (GPS) and total station) can provide slope stability information from point deformation measurements; however, their sparse spatial distributions reduce their applicability in large-scale areas with rugged terrain. The satellite-based interferometric synthetic aperture radar (InSAR) has proven to be a powerful tool to monitor and analyze catastrophic slope failure in vast and inaccessible regions. The microwave measurement enables InSAR to monitor ground movement in centimeter level with high spatial resolution and broad coverage under all weather conditions [9,10]. Benefiting from the extensive archive of SAR acquisitions covering the same area, multi-temporal InSAR (MTInSAR) have been developed in recent years to overcome the limitations associated with conventional InSAR techniques by exploiting time-series information of SAR acquisitions, all of which basically fall into two categories, namely persistent scatterer interferometry (PSI) [11][12][13][14] and short baseline subset (SBAS) [15][16][17][18]. PSI exploits phase information from point-like persistent scatterers (PSs) that keep high coherence in long observation period and usually correspond to man-made objects and large rock outcrops. This method uses N + 1 single-look SAR images to form N single-master interferograms with good performance in urban or peri-urban areas [19,20]. SBAS, on the other hand, extracts phase information from distributed scatterers (DSs) that have moderate coherence in interferograms with short spatial and temporal baselines. The DSs are referred to as targets that contain a sum of random scatterers which are spatially multilooked to improve the signal-to-ratio (SNR) during SBAS processing. The time-series InSAR processing has been widely applied in mapping ground deformation associated with seismic activities [21], volcanism [22], glacier movement [23], underground resource exploitation [24], and landslides [5,25,26], etc. However, InSAR is less effective in low coherence environments. Areas covered by dense vegetation are usually excluded from the analysis, resulting in gaps in an InSAR derived deformation map.
Despite the sustainable development of InSAR algorithm, InSAR technique is still less effective in low-coherence environments [27,28]. Areas covered by dense vegetation are usually excluded from the analysis, resulting in gaps in an InSAR derived deformation map. Landslides have severe impacts on local residents in Taiwan due to the unique geological conditions and high annual rainfall [29,30]. Highly vegetated areas and rugged terrain in Taiwan pose great challenge to conventional InSAR techniques for quantitative hazard and risk assessment.
In this study, we modify the existing MTInSAR algorithm to process phase information from both PSs and DSs and thereby improve the spatial coverage in environments with low coherence. Phase denoising and phase unwrapping are treated as an integrated process to avoid extra estimation error due to necessary approximations and presumptions in each individual step. Systematic errors are corrected step-by-step, so improve the final deformation precision. We apply the method to characterize the time series slope displacements prior to the landslide failure of Freeway No. 3 in northern Taiwan. In addition, we address whether the long-term slope movement is correlated with precipitation. This study is expected to demonstrate the applicability of InSAR technique for safety inspections of geotechnical structures in mountainous areas.  (Figure 1). The slope originally ranged between 12 • and 20 • with a northwest to southeast aspect at an elevation of 110 m to 162 m. The escarpment forms a cuesta landscape, in which one slope is gentle and follows the dip of the bedding (the sliding surface), and the other slope is steep and rugged due to the differential erosion of the sandstone and shale interbeds.

Study Area and Datasets
Remote Sens. 2020, 12, 2403 3 of 19 The exposed strata are from the Miocene (Tertiary) Taliao Formation, which is composed of massive sandstone, interbeds of sandstone and shale, and thick-bedded sandstone intercalated with shale ( Figure 2). The strike orientations are approximately from N20 • E/15 • SE to N40 • E/15 • SE, and that brittle structure indicates as a conjugate shear fractures aligned N66 • W/68 • NE and N45 • W/73 • SW. The slope on the right side of the southbound lane is a dip slope. The boring core sample from the dip slope shows that the slope's 10 m-thick surface layer is thick-bedded sandstone occasionally intercalated with shale and sits atop interbeds of sandstone and shale. In addition, the average rock quality designation (RQD) at the boundary of the sandstone and the shale is no more than 25%, implying that the rock quality of this section is poor [3].
Remote Sens. 2020, 12, x FOR PEER REVIEW 3 of 19 landscape, in which one slope is gentle and follows the dip of the bedding (the sliding surface), and the other slope is steep and rugged due to the differential erosion of the sandstone and shale interbeds. The exposed strata are from the Miocene (Tertiary) Taliao Formation, which is composed of massive sandstone, interbeds of sandstone and shale, and thick-bedded sandstone intercalated with shale ( Figure 2). The strike orientations are approximately from N20°E/15°SE to N40°E/15°SE, and that brittle structure indicates as a conjugate shear fractures aligned N66°W/68°NE and N45°W/73°SW. The slope on the right side of the southbound lane is a dip slope. The boring core sample from the dip slope shows that the slope's 10 m-thick surface layer is thick-bedded sandstone occasionally intercalated with shale and sits atop interbeds of sandstone and shale. In addition, the average rock quality designation (RQD) at the boundary of the sandstone and the shale is no more than 25%, implying that the rock quality of this section is poor [3].  Known as a freeway excavation site, this area has an average excavation depth of 17 m as its original ground elevation was 120 m, while the designed road elevation is 103 m. Together with the excavation of side slopes, the site's total excavation height is approximately 35 m, resulting in the exposure of not only thick-bedded sandstone, but also shale intercalated with sandstone and shale and sandstone interbeds. The collapsed slope is located within the freeway excavation site and was originally excavated at 30 to 45 degrees; The upper section was a natural slope excavated at 30 degrees, whereas the lower section was a two-stage 10-m slope excavated at 45 degrees with failure protection in the form of prestressed rock anchors and precast reinforced concrete (RC) grid beams. Known as a freeway excavation site, this area has an average excavation depth of 17 m as its original ground elevation was 120 m, while the designed road elevation is 103 m. Together with the excavation of side slopes, the site's total excavation height is approximately 35 m, resulting in the exposure of not only thick-bedded sandstone, but also shale intercalated with sandstone and shale and sandstone interbeds. The collapsed slope is located within the freeway excavation site and was originally excavated at 30 to 45 degrees; The upper section was a natural slope excavated at 30 degrees, whereas the lower section was a two-stage 10-m slope excavated at 45 degrees with failure protection in the form of prestressed rock anchors and precast reinforced concrete (RC) grid beams.
Influenced by the subtropical-monsoon climate, the study area suffers typhoons in the summer and the northeast monsoon in the winter. On average, there are 190 rainy days per year, yielding an annual precipitation of more than 2000 mm [29,30]. Most of the rainfall occurs in the typhoon season, Remote Sens. 2020, 12, 2403 4 of 19 which is from July to September. Due to mountainous terrain with steep slopes and fractured bedrock, the rainfall associated with these tropical storms increases the risk of landslides and debris flows. Influenced by the subtropical-monsoon climate, the study area suffers typhoons in the summer and the northeast monsoon in the winter. On average, there are 190 rainy days per year, yielding an annual precipitation of more than 2000 mm [29,30]. Most of the rainfall occurs in the typhoon season, which is from July to September. Due to mountainous terrain with steep slopes and fractured bedrock, the rainfall associated with these tropical storms increases the risk of landslides and debris flows.

The Landslide
The slope failure of National Freeway No. 3 landslide is a typical example of translational dipslope slide as the slump moved eastwards along the shale bedding plane in the sandstone and shale interbeds [4,31]. No rainfall or earthquake was recorded in weeks prior to the event. Heavily covered in vegetation (shown in Figure 3), the slide mass mainly consists of the surface soil layer and the broken weathered thick-bedded sandstone occasionally intercalated with shale. The 200,000 m 3 of slide mass not only damaged the Dapu Bridge overpass, but also interrupted traffic and killed 4 people as its 200 m-long and 150 m-wide slump covered the six-lane freeway.

The Landslide
The slope failure of National Freeway No. 3 landslide is a typical example of translational dip-slope slide as the slump moved eastwards along the shale bedding plane in the sandstone and shale interbeds [4,31]. No rainfall or earthquake was recorded in weeks prior to the event. Heavily covered in vegetation (shown in Figure 3

Satellite SAR Data
An archived SAR dataset of L-band ALOS/PALSAR data is used to retrieve the surface movement before the disaster. Compared with other SAR sensors with shorter radar wavelength, the L-band PALSAR data have better performance in vegetated environments. A total of 15 PALSAR images were acquired from an ascending track during the period from January 2007 to March 2010 with a minimum revisiting cycle of 46 days. Figure 4 illustrates the combination of 105 interferometric pairs from 15 SAR images. Considering the narrow region of the collapsed slope, full resolution interferograms are generated with a spatial resolution of 10 m. In addition, a 1-arc-second ALOS World 3D 30 m (AW3D30m) digital surface model (DSM) [32] is used to remove the topographic phase contribution to the interferometric phase.

Satellite SAR Data
An archived SAR dataset of L-band ALOS/PALSAR data is used to retrieve the surface movement before the disaster. Compared with other SAR sensors with shorter radar wavelength, the L-band PALSAR data have better performance in vegetated environments. A total of 15 PALSAR images were acquired from an ascending track during the period from January 2007 to March 2010 with a minimum revisiting cycle of 46 days. Figure 4 illustrates the combination of 105 interferometric pairs from 15 SAR images. Considering the narrow region of the collapsed slope, full resolution interferograms are generated with a spatial resolution of 10 m. In addition, a 1-arc-second ALOS World 3D 30 m (AW3D30m) digital surface model (DSM) [32] is used to remove the topographic phase contribution to the interferometric phase. Remote Sens. 2020, 12, x FOR PEER REVIEW 6 of 19

Methodology
To improve the density and reliability of deformation signal in densely vegetated environments, we modify the existing MTInSAR analysis framework to retrieve the spatial and temporal evolution of landslide movement. The modified framework consists of three steps, i.e., coherent point selection, interval phase estimation, and systematic error correction. Figure 5 displays the flowchart of the MTInSAR analysis. Particularly, coherent phase observations are extracted from the combination of selected PS and DS measurements. The interval phase estimation is performed at arcs then converted to points. Then, the interval phase maps are corrected by removing topographic residuals and atmospheric delays. Finally, the time-series displacements are derived from the conversion of the remaining phase component. In the following the detailed procedure of the modified MTInSAR processing will be explained.

Methodology
To improve the density and reliability of deformation signal in densely vegetated environments, we modify the existing MTInSAR analysis framework to retrieve the spatial and temporal evolution of landslide movement. The modified framework consists of three steps, i.e., coherent point selection, interval phase estimation, and systematic error correction. Figure 5 displays the flowchart of the MTInSAR analysis. Particularly, coherent phase observations are extracted from the combination of selected PS and DS measurements. The interval phase estimation is performed at arcs then converted to points. Then, the interval phase maps are corrected by removing topographic residuals and atmospheric delays. Finally, the time-series displacements are derived from the conversion of the remaining phase component. In the following the detailed procedure of the modified MTInSAR processing will be explained.

Methodology
To improve the density and reliability of deformation signal in densely vegetated environments, we modify the existing MTInSAR analysis framework to retrieve the spatial and temporal evolution of landslide movement. The modified framework consists of three steps, i.e., coherent point selection, interval phase estimation, and systematic error correction. Figure 5 displays the flowchart of the MTInSAR analysis. Particularly, coherent phase observations are extracted from the combination of selected PS and DS measurements. The interval phase estimation is performed at arcs then converted to points. Then, the interval phase maps are corrected by removing topographic residuals and atmospheric delays. Finally, the time-series displacements are derived from the conversion of the remaining phase component. In the following the detailed procedure of the modified MTInSAR processing will be explained.

Coherent Point Selection
It is desired to retrieve ground coherent points as many as possible for the purpose of sufficiently extracting slow slip displacement over vegetated area. To this end, the basic observations in this method contain both PS and DS measurements. For PS, its phase stability can be approximated by temporal amplitude stability [11], which is expressed by amplitude dispersion index (ADI), where µ A and σ A represent the mean and standard deviation of time-series amplitudes for a given pixel, respectively. Pixels with ADI less than 0.4 are initially considered as PS candidates [11,12]. In contrast to PS, DS maintains moderate coherence during one or several intervals of SAR acquisitions. For a given pixel within an interferometric pair, its complex coherence value is calculated as, where s m (t) and s n (t) represent the complex phase values of SAR image pair, * stands for complex conjugate, the coherence magnitude γ m,n defines the degree of coherence, and φ m,n represents the spatially averaged interferometric phase, L is the number of statistically homogeneous pixels (SHPs) surrounding the central pixel. The SHPs for the central pixel can be identified through statistical similarity test (e.g., Kolmogorov-Smirnov (KS) and Anderson-Darlington (AD) test) based on time-series amplitudes [33][34][35]. Points with average temporal coherence larger than 0.25 and number of SHPs larger than 20 are initially selected as DS candidates [36]. The original phases of selected DS points are substituted with the spatially averaged phases. To preserve the phases of point-wise PSs, the DS points that coincide with PSs points are dropped. Noted that, the coherence of PS is always equal to one in all interferograms [33]. Then the selected PSs and DSs are combined for joint process in the following analysis. Their phase and coherence values serve as the basis of the functional and stochastic models of the interval phase estimation, respectively.

Interval Phase Estimation
Once the coherent points have been selected, phase differences at arcs constructed from Delaunay triangulation form the observations for interval phase estimation. Since the arc-based observations between two neighboring coherent points can largely reduce 2π ambiguities in phase measurements [18,37], we first estimate the phase intervals at arcs. Assuming there are N + 1 SAR images and M interferometric pairs, the linear system for arcs without considering phase ambiguities is expressed as [38] where ∆φ is a vector containing the phase differences between two adjacent pixels from M interferometric combinations, ∆ϕ is the estimated phase difference at arc for N SAR acquisition intervals, A is the design matrix connecting SAR acquisition interval and interferogram, R is a stochastic vector with an expectation of zero. P is equivalent weight matrix, which adjusts the influence of observations to fit the actual accuracy of the corresponding observations [39]. It is expressed as the combination of coherence fisher information matrix (FIM) [40] and down-weighting functions [41]. We solve the linear system (3) by iteratively updating weights according to the residuals until the model is fitted, so the effect from outliers (i.e., phase ambiguities) is suppressed [42]. In case the percentage of outliers exceeds of the upper limit of model and results in corrupted estimation, the arc measurements Remote Sens. 2020, 12, 2403 8 of 19 need to be removed to avoid contaminating subsequent estimation. Such arcs can be identified by equivalent temporal coherence (ETC), where tr(·) represents the operator to calculate trace of matrix. Arcs with γ P less than a given threshold (e.g., 0.85) are considered low-quality and are discarded [33,43]. We perform the estimation of (3) and (4) on an arc-by-arc basis until the phase chain at all arcs are obtained. Isolated points are detected and removed to avoid rank deficiency in design matrix A. Then we convert the estimated phase differences at arcs to unwrapped phases at points by spatially weighted integration, where C is network matrix relating points and arcs, which is composed by the combination −1, 0 and 1. W is the weight matric of arcs, of which diagonal element is determined by the posterior variance from Equation (3). We refer the readers who are interested in the details of unwrapped interval phase retrieval to Liang, et al. [38].

Systematic Error Correction
After obtaining interval phases at coherent points, the consecutive phase map can be thought as a linear mixture of statistically independent signals, δϕ = δϕ defo + δϕ topo + δϕ atm + δϕ orb (6) where δϕ de f o represents the phase component due to ground movement, δϕ topo represents the phase residual caused by digital elevation model (DEM) inaccuracy, δϕ atm is the phase screen due to atmospheric delays, δϕ orb is the phase ramp induced by orbital inaccuracy. Because phase ambiguities have been removed, we only need to separate the error components from the interval phase maps, thus the remaining phase component can be easily converted to displacement. Firstly, the DEM error can be corrected through the linear relation between spatial baseline and phase residual. Since the linear estimation is prone to the discrepancy between the modeled and real deformation time series [44], we employ a component-based method to correct the topographic error. The detailed analysis and test of the method has been discussed in Liang, et al. [45]. The basic idea stems from the fact that the topographic error has a deterministic spatial pattern which is distinct with other signals and merely correlated with spatial baselines. By using independent component analysis (ICA), the spatial pattern due to topographic error can be identified without requiring a priori deformation model.
After correcting the DEM error, there are deformation, orbital error and atmospheric delays remained in the interval phases. The height-dependent atmospheric delays can be corrected by a quadtree aided joint model, which has been discussed and tested with simulated and real data in Liang, et al. [46]. We recall the method and extend it to simultaneously account for height-dependent atmospheric delays and orbital error. For a given SAR acquisition, the phase ramp due to height-dependent atmospheric delays and orbital error with respect to a reference image is expressed as where x and y are the rang and azimuth coordinate, respectively. h represents the elevation for SAR pixel, k 1 to k 4 are the coefficients which correlate orbital error and height-dependent atmospheric delay terms to the phase ramps. φ 0 is the intercept with respect to the reference point, which can be neglected where Φ atm+orb represents (J × N) × 1 vector including phase components caused by height-dependent atmospheric delays and orbital error, P atm+orb contains all orbital error and height-dependent atmospheric delay coefficients with respect to the reference SAR acquisition and reference point, B atm+orb is the design matrix connecting SAR intervals and acquisitions for the unknown parameters. Apart from the height-dependent atmospheric delays and orbital error, we model the displacement as where λ is the radar wavelength, T is the temporal baseline of interval phase, v is the deformation rate. There are J + 1 points and N + 1 SAR acquisitions, the phase component contributed from the modelled deformation is expressed as [47], where Φ de f o represents (J × N) × 1 vector including phase component caused by modelled displacements, P de f o contains the parameters of deformation rate of all points with respect to the reference point. B de f o is the design matrix connecting SAR intervals and acquisitions for the deformation parameters. Combining (8) and (10), we derive the joint model as where Φ represents (J × N) × 1 vector including all phase components in the joint model, R is the phase residual vector. By solving (11), the parameters of height-dependent atmospheric delays, orbital error and deformation rate can be obtained simultaneously. We further subtract the height-dependent atmospheric delays and orbital error from the interval phase maps, so the time-series displacements are directly transformed from the remaining phase component.

Phase Reconstruction and Systematic Error Correction
The improved MTInSAR approach is applied to the ALOS/PALSAR data as detailed in Section 2.2. We use ADI to select 6965 PS points and use mean coherent map to select 82,855 DS points from 105 single-look interferograms. The phase measurements from the selected PS and DS points are combined as coherent observations in interval phase estimation. After rejecting low-quality arcs and removing isolated points, we obtain 14 interval phase maps with 80,422 points. To illustrate the performance of the proposed method on noise reduction, Figure 6 displays four examples of single-look consecutive interferograms before and after phase interval estimation and systematic error correction. From the first row in Figure 6 we can see that the phases in the original single-look interferograms are deteriorated from decorrelation noise. Although the L-band SAR sensor has comparable performance at maintaining coherence in rural area, the temporal decorrelation due to the unfavorable geology and dense vegetation make the interferometric phase noisy and hampers the retrieval of landslide movement. The second and third rows in Figure 6 display the consecutive interferograms in wrapped and unwrapped form after phase estimation from the combined measurements of PSs and DSs. It is clear that the decorrelation noise is effectively filtered, yielding a more coherent phase map for deformation retrieval. The interval phase maps are further processed by correcting topographic error and height-dependent atmospheric delays and orbital error. The comparison between the third and fourth rows in Figure 6 indicates that, besides decorrelation noise, the phases in the interferograms also have contributions from the DEM errors and phase ramps. By subtracting the estimated error phase components, the unexpected phase variations are largely reduced, indicating the effectiveness of systematic error correction in the proposed MTInSAR processing.
After performing interval phase estimation and error mitigation, we convert the remaining phase component into a displacement time series and calculate the deformation rate. The average point density reaches approximately 3251 points/km 2 . The high-density deformation measurements over the study area help to better identify and analyze the slow-moving landside before the final slope failure.  Figure 7 shows the annual deformation rate map of the investigated slope area during the 2007-2010 period. Although the studied slope was covered by dense vegetation (Figure 3), the results show considerable deformation distributed along the slope that subsequently collapsed. The deformation is concentrated on the upper part of the slope, with the average LOS velocity reaching −30mm/year over the investigated period. Given the ascending ALOS imaging geometry and the slope facing an approximately east-southeast direction, the negative LOS displacement i.e., moving away from the  Figure 7 shows the annual deformation rate map of the investigated slope area during the 2007-2010 period. Although the studied slope was covered by dense vegetation (Figure 3), the results show considerable deformation distributed along the slope that subsequently collapsed. The deformation is concentrated on the upper part of the slope, with the average LOS velocity reaching −30mm/year over the investigated period. Given the ascending ALOS imaging geometry and the slope facing an approximately east-southeast direction, the negative LOS displacement i.e., moving away from the satellite, indicates downslope movement. The distribution of downslope movement is consistent with the region that collapsed afterward, which is enclosed by the solid red triangle. Figure 8 shows the time-series LOS displacement maps of the studied slope. It can be observed that there are two sliding acceleration periods, i.e., July 2007 to January 2008, and July 2009 to March 2010, which can possibly be ascribed to the precursor of the slope failure that occurred on 25 April, 2010. Specifically, the larger amplitude downslope movement occurs on the upper part of the slope, and the amplitude is positively correlated with the slope elevation. We also note that, the area on the left side shows positive deformation values. The uplift movement can be explained from two aspects. First, the orientation of the left slope is different from that of the central slope. Or, alternatively, the uplift is caused by the horizontal movement since the InSAR measurements are in the LOS direction.

Spatial Pattern and Temporal Behavior of Slope Movement
Remote Sens. 2020, 12, x FOR PEER REVIEW 11 of 19 positively correlated with the slope elevation. We also note that, the area on the left side shows positive deformation values. The uplift movement can be explained from two aspects. First, the orientation of the left slope is different from that of the central slope. Or, alternatively, the uplift is caused by the horizontal movement since the InSAR measurements are in the LOS direction. Figure 9 shows the time series displacements of four selected points in Figure 7. Point 1 and Point 2, located in the upper section with a relatively steep slope, move downwards with mean velocities of −18.6mm/year and −21.3mm/year, respectively. Point 3, located in the middle section of the slope, has a mean downslope velocity of −15.0mm/year. While Point 4 has an overall mean velocity of −0.1mm/year, there is an apparent seasonal oscillation in the time series, possibly as the result of interaction between the slope and the freeway. All the four points show complex temporal behaviors such that the movement cannot be simply scaled through time during the period of observation, indicating the landslide movement was possibly controlled by several factors. Active ground surface deformation identified in the InSAR data can be classified into three types according to ground features before the collapse: first, the primitive forest located on the upper slope, where Point 1 (crown of the slope with obvious tension cracks) and 2 (center of the sliding block) are located; second, the artificial slope, where a hairpin road can be clearly identified in Figure 7; third, the boundary of the primitive forest and artificial slope, where Point 4 (rock anchoring area) is located.  location of selected points in (Figure 9). The black dash line marks the location displacement profile shown in Figures 2 and 10.   Figure 9 shows the time series displacements of four selected points in Figure 7. Point 1 and Point 2, located in the upper section with a relatively steep slope, move downwards with mean velocities of −18.6mm/year and −21.3mm/year, respectively. Point 3, located in the middle section of the slope, has a mean downslope velocity of −15.0mm/year. While Point 4 has an overall mean velocity of −0.1mm/year, there is an apparent seasonal oscillation in the time series, possibly as the result of interaction between the slope and the freeway. All the four points show complex temporal behaviors such that the movement cannot be simply scaled through time during the period of observation, indicating the landslide movement was possibly controlled by several factors. Active ground surface deformation identified in the InSAR data can be classified into three types according to ground features before the collapse: first, the primitive forest located on the upper slope, where Point 1 (crown of the slope with obvious tension cracks) and 2 (center of the sliding block) are located; second, the artificial slope, where a hairpin road can be clearly identified in Figure 7; third, the boundary of the primitive forest and artificial slope, where Point 4 (rock anchoring area) is located.

Interpretation of the Freeway No. 3 Landslide Movement
To further visualize the spatial and temporal pattern of the landslide, we extracted the profile from the average deformation rate and the cumulative displacements in the sliding zone. The profile, which is denoted by black dashed line in Figure 7, has a head-to-toe length of approximately 480 m and elevation change of approximately 70 m in the northeast-southwest direction. From Figure 10, we can see that the downslope movement decreases along the slope and even uplift near the road, indicating that the retaining wall of the side slope played a role to restraining the sliding mass. The area with active movement and the highest LOS velocity is located on the upper edge of the side slope, which is between the head scarp at the hilltop and Freeway No. 3 at the base of the slope.
According to the on-site measurement at several different positions, the average strike of the bedding plane is approximately from N10 • E to N32 • E with a dip angle between 12 • SE to 15 • SE. On the other hand, the strike of both slopes to either side of the landslide is N33 • E, and its bedding plane dip angle is only approximately 10 • . This finding not only means that the collapse mode involves a detached block slide along weak planes wherein the dip direction of a weak plane is nearly parallel to the slope surface. It also implies that the failure mechanism was triggered by a translational slide [3,4,31].
The field survey results reveal that the dip slope failure surface was along a shale bedding plane within sandstone and shale interbeds (Figure 2). To the south, the sliding block failed along the upper edge of the first layer of the anchor-reinforced slope, while in the north, the slide began near the base of the first layer of the anchor-reinforced slope. Although exposed slickensides were apparent, no water was extruded from the sliding block or from the edge of sliding zone, and the steel strands from rock anchors were found locally within the region.
Upon examination of the broken rock anchor heads that slid to the top of the hill on the eastern side of the slope, it was found that some steel strands remaining inside the anchor heads were seriously corroded and torn. Because the fractures in these steel strands are mostly corroded and tainted, they are considered old fractures. It is therefore inferred that the said steel strands may have broken a significant amount of time before this incident. The slow downslope movement could be responsible for the gradual process of anchorage loss before the slope collapsed.
we can see that the downslope movement decreases along the slope and even uplift near the road, indicating that the retaining wall of the side slope played a role to restraining the sliding mass. The area with active movement and the highest LOS velocity is located on the upper edge of the side slope, which is between the head scarp at the hilltop and Freeway No. 3 at the base of the slope.
According to the on-site measurement at several different positions, the average strike of the bedding plane is approximately from N10°E to N32°E with a dip angle between 12°SE to 15°SE. On the other hand, the strike of both slopes to either side of the landslide is N33°E, and its bedding plane dip angle is only approximately 10°. This finding not only means that the collapse mode involves a detached block slide along weak planes wherein the dip direction of a weak plane is nearly parallel to the slope surface. It also implies that the failure mechanism was triggered by a translational slide [3,4,31].

Correlation with Rainfall
It has been well-documented that the ground movement of landslide evolves in response to variation of precipitation [5,6,48]. Infiltration of surface water causes an increase in pore-water pressure and the driving force of slope failure, which in turn decreases the effective normal stress and resisting friction along the shear zone [49,50]. To investigate the temporal correlation between landslide movement and precipitation in this study, we compare precipitation and the displacement time series obtained from our InSAR measurements. We use precipitation data from five surrounding meteorological stations within 15 km of the landslide location and calculate their mean value to use as the estimated rainfall on the studied slope. Then, we calculate the average amount of incremental displacement for each time period from the enclosed landslide region in Figure 7. The correlation between the incremental displacement from InSAR and the precipitation measurements ( Figure 11) shows that the landslide has intermittent movement in the time series that correlates with rainfall. It is observed that the study area generally has high rainfall in October and December and the active movement clearly accelerates during these periods. Although the correlation is not perfect, a similar temporal pattern can be found between precipitation and incremental displacement, especially in the period from January 2007 to October 2008. According to the InSAR time-series analysis results, the landslide activity of the slope can be divided into four stages. The first stage is between January and July 2007 with insignificant downward movement. Then, between July 2007 and January 2008, the slope displacement rate obviously increased with a displacement of -15mm. Between January 2008 and July 2008, the slope became stabilized again with insignificant displacement. As ALOS did not capture any image of Taiwan between October 2008 and July 2009, it is difficult to evaluate the slope deformation during this period. Nevertheless, the slope's overall displacement rate still reached -10 mm between July 2008 and July 2009. Finally, the last stage was between July 2009 and March 2010, where the slope's overall displacement rate was -20 mm. Comparing with the previous two stages (January to July 2007, and January to July 2007), the overall slope deformation continued in October 2009 during the dry season. The underlying mechanism can be summarized as follows. The permeable and weathered upper layer of sandstone provides favorable condition for the infiltration of rainfall into the lower shear layer, which is relatively impermeable and hampers drainage water. The accumulation of water at the stratigraphic boundary reduces the resisting friction and causes the corrosion of ground anchors. When the force of the sliding rock mass exceeds the supporting force from the rock anchors and resisting friction, the slope is inclined to move downwards, triggering a major landslide.
Remote Sens. 2020, 12, x FOR PEER REVIEW 15 of 19 Figure 11. Comparison of the incremental downslope displacement (blue lines) with the average monthly precipitation (gray histogram).

Performance of InSAR Technique on Detecting and Monitoring Landslides
In this study, we use 15 ALOS/PALSAR images to generate full combination of 105 interferograms, which provide sufficient observations for the deformation analysis. The following three factors contribute to the improvement of the density and reliability of deformation signal in densely vegetated environments.


Coherent point selection: although PS has stable scatterering behavior over long-term periods, the sparse distribution of PS in rural areas compromises the efficiency of landslide displacement retrieval. To increase the density of coherent measurements, not only PS points are included in the processing analysis, but also DS points, which are widely distributed in mountainous regions and correspond to ground targets with high sensitivity to spatial and temporal baselines.  Interval phase estimation: the objective of this procedure is to obtain sufficient optimized phase measurements, which involves phase denoising and phase unwrapping. In traditional processing flow, phase denoising and phase unwrapping are carried out separately (e.g., phase denoising is achieved by phase triangulation [33], while phase unwrapping is implemented by minimum cost flow (MCF) [51]). However, the separated processes will introduce extra estimation error due to necessary approximations and presumptions in each step [52]. On the contrary, in this study, we apply an arc-based phase estimation approach to jointly conduct phase denoising and phase unwrapping. The linear formulation via this method ensures the estimation efficiency and quality of the phase optimization.

Performance of InSAR Technique on Detecting and Monitoring Landslides
In this study, we use 15 ALOS/PALSAR images to generate full combination of 105 interferograms, which provide sufficient observations for the deformation analysis. The following three factors contribute to the improvement of the density and reliability of deformation signal in densely vegetated environments.
• Coherent point selection: although PS has stable scatterering behavior over long-term periods, the sparse distribution of PS in rural areas compromises the efficiency of landslide displacement retrieval. To increase the density of coherent measurements, not only PS points are included in the processing analysis, but also DS points, which are widely distributed in mountainous regions and correspond to ground targets with high sensitivity to spatial and temporal baselines.

•
Interval phase estimation: the objective of this procedure is to obtain sufficient optimized phase measurements, which involves phase denoising and phase unwrapping. In traditional processing flow, phase denoising and phase unwrapping are carried out separately (e.g., phase denoising is achieved by phase triangulation [33], while phase unwrapping is implemented by minimum cost flow (MCF) [51]). However, the separated processes will introduce extra estimation error due to necessary approximations and presumptions in each step [52]. On the contrary, in this study, we apply an arc-based phase estimation approach to jointly conduct phase denoising and phase unwrapping. The linear formulation via this method ensures the estimation efficiency and quality of the phase optimization. • Systematic error correction: Multiple error sources exist in phase measurements despite the implementation of phase optimization. In this study, topographic residual and height-dependent atmospheric delays are the most prominent ( Figure 6). The component-based method provides an efficient way to correct the topographic residual without requiring a priori deformation model. Although the joint model is time-consuming, it simultaneously estimates height-dependent atmospheric delays, orbital error and deformation signal, avoiding estimation bias caused by confounding processes [47].
The developed MTInSAR method demonstrates the applicability of InSAR technique in determining landslide boundary and assessing slope instability, especially where ground truth data are scarce. However, some limitations of this method still exist. Firstly, although L-band ALOS\PALSAR data has good performance in the presence of dense vegetation and phase optimization further reduces decoherence noise, it cannot be extended to all types of ground surfaces and SAR sensors, especially for higher frequency band data and areas with rapid surface changes which exceeds the maximum detectable phase range (i.e., 1/4 of the radar wavelength) [27]. Secondly, single track InSAR measurements can only retrieve the movement component along the LOS direction, and it shows low sensitivity to the movement that is perpendicular to the LOS direction. This could be done with the aid of SAR amplitude information [53] or corner reflectors that provide high electromagnetic reflectivity over vegetated slopes. Meanwhile, a higher temporal sampling could increase the chance of capturing more detailed deformation variation and finding the precursor of the failure process. Secondly, single track InSAR measurements can only retrieve the movement component along the LOS direction, and it shows low sensitivity to the movement that is perpendicular to the LOS direction.
To remedy the lack of information, exploiting additional assumptions on ground movement or combing ascending and descending tracks could potentially recover the actual slope displacement.

Conclusions
An optimized MTInSAR approach based on the integrated use of PSs and DSs has been developed in this study. The approach jointly deals with phase denoising and phase unwrapping, so avoid extra estimation error in each individual step. Systematic errors are corrected step-by-step according to their spatial and temporal characteristics. We apply this approach to analyze the precursor movement of the landslide failure along Freeway No. 3 in northern Taiwan. The analysis results indicate that before the catastrophic landslide failure, variable downslope movement occurred along this dip slope. The spatial and temporal variations in deformation are in good agreement with previous studies from field and laboratory investigations. The temporal behaviors further reveal that the time-series displacement of the slope is correlated with the local precipitation. This study demonstrates the informative and beneficial contribution provided by InSAR technique in landslide research, especially where ground surface is covered by dense vegetation and ground truth data are scarce. One insight from this slope failure analysis is that landslide failure, either along natural or artificial slopes, have long-term signs of subtle movement that can be detected by satellite radar several months before the final collapse. Continuous safety inspection and monitoring using InSAR techniques could be considered as an efficient way to examine the life cycles of such geotechnical structures.