Mapping Two-Dimensional Deformation Field Time-Series of Large Slope by Coupling DInSAR-SBAS with MAI-SBAS

Mapping deformation field time-series, including vertical and horizontal motions, is vital for landslide monitoring and slope safety assessment. However, the conventional differential synthetic aperture radar interferometry (DInSAR) technique can only detect the displacement component in the satellite-to-ground direction, i.e., line-of-sight (LOS) direction displacement. To overcome this constraint, a new method was developed to obtain the displacement field time series of a slope by coupling DInSAR based small baseline subset approach (DInSAR-SBAS) with multiple-aperture InSAR (MAI) based small baseline subset approach (MAI-SBAS). This novel method has been applied to a set of 11 observations from the phased array type L-band synthetic aperture radar (PALSAR) sensor onboard the advanced land observing satellite (ALOS), spanning from 2007 to 2011, of two large-scale north–south slopes of the largest Asian open-pit mine in the Northeast of China. The retrieved displacement time series showed that the proposed method can detect and measure the large displacements that occurred along the north–south direction, and the gradually changing two-dimensional displacement fields. Moreover, we verified this new method by comparing the displacement results to global positioning system (GPS) measurements.


Introduction
Landslide is one of the widely distributed and frequently occurring geological disasters, which commonly causes many casualties as well as huge property loss [1,2]. Monitoring the temporal and spatial variations of slopes, such as large-scale hill slopes and engineering slopes, is essential for dynamic analysis, early-warning, and risk assessment [3][4][5]. Satellite radar remote sensing system has more significant advantages in the precise analysis of landslides over large areas than ground-based monitoring systems, such as digital leveling instruments, total stations, global navigation satellite system (GNSS) receivers, extensometers, inclinometers, and so on. Therein, the spaceborne synthetic aperture radar (SAR) data acquisition and processing system, being an advanced radar remote sensing system, could provide valuable temporal and spatial evolution information in all weather conditions for landslide monitoring [6][7][8].
A typical characteristic of landslide is that its displacements occur in both vertical and horizontal directions. However, the conventional differential SAR interferometry (DInSAR) technique can only measure one-dimensional (1-D) displacement along the line-of-sight (LOS) direction [9,10]. It means that any component of motion orthogonal to the LOS cannot be detected by DInSAR approach. For the nearly north-south flight pattern of current commercial SAR satellites, DInSAR is only sensitive to the displacement along the up-down and east-west directions [11,12]. Therefore, it is almost unable to measure the displacement occurring in the north-south direction with the conventional DInSAR method.
To cope with the problem of displacement monitoring along the north-south direction, Michel et al. [13] proposed an offset-tracking method utilizing radar amplitude information to measure azimuthal offsets. Offsets reflect the difference in position in the two radar images of a given point of the ground. This difference can be converted to ground displacements. However, this method involves a sub-pixel correlation of the radar images and is computationally expensive [14]. Furthermore, the accuracy of the offset-tracking method is sensitive to the spatial resolution of the SAR data and the precision of the co-registration [15].
In 2006, a significant improvement in measuring along-track displacement was made with a new approach: multiple-aperture InSAR (MAI) [16]. The along-track component of displacement can be measured by using sub-aperture processing. The phase difference between the forward-and backward-looking interferograms is proportional to any along-track displacement. The MAI method is more accurate than the offset-tracking approach with less computation time, and it is easy to implement using a conventional InSAR software. Subsequently, Gourmelen et al. utilized the MAI method to measure the ice surface velocity on two ice caps in Iceland [17]. Recently, the MAI method was used to estimate the three-dimensional (3-D) displacement caused by earthquakes and volcanoes by combining multiple-aperture interferometry and conventional interferometry [18,19]. Moreover, the stacking MAI technique has been used to measure the slow-moving azimuthal displacements [20].
Compared to DInSAR, MAI can be generally used in areas with large-scale displacements because it more sensitive to phase errors [21]. This limitation is due to the fundamental tradeoff between sensitivity and signal to noise ratio (SNR) in partial aperture processing. Therefore, whether MAI can be used to monitor landslide is still a question. In this paper, we try to perform a feasibility check for large-scale landslide monitoring by combining DInSAR and MAI. Firstly, we intend to evaluate the capability of the MAI technique to map the horizontal displacement of two large-scale slopes with L-band SAR data. Secondly, we explore the possibility of extending MAI interferograms into time series by introducing the main idea of small baseline subset (SBAS) method. Thirdly, we intend to combine MAI-SBAS and conventional DInSAR-SBAS to retrieve the multi-temporal displacement field of two large-scale landslides. To reach the aforementioned goals, 11 acquisitions of L-band SAR data of Fushun west open-pit mine (FWOPM), which is the largest open-pit in Asia, from January 2007 to January 2011, were used to reveal the displacement field of landslides caused by coal mining from the overall perspective via radar remote sensing technique.
This paper is organized as follow. Firstly, the study area, including its social economic environment and landslide background, as well as the available SAR data, is introduced. Secondly, the methodology used in this study is described, highlighting the algorithm procedures of DInSAR-SBAS and MAI-SBAS and the combination approach for both results. Subsequently, we present the results obtained with the DInSAR-SBAS method, the MAI-SBAS method and the combination method, respectively, applied in the study area. The following section is focused on the comparison between the achieved InSAR results and the global positioning system (GPS) measurements, and we further compared our results with the high spatial resolution optical image and other studies. The last section is dedicated to the conclusions, summarizing the main findings of this study.

Study Area and SAR Data
Located in Fushun city, Liaoning province, Northeast of China, FWOPM has a long history of more than one hundred years and makes important contributions to China's economic development. After a hundred years mining, FWOPM has produced the largest open-pit mine in Asia (Figure 1), with a length of 6600 m east-to-west, width of 2200 m north-to-south, and an average depth approximately 420-520 m, occupying the lowest terrestrial point in mainland China [22,23]. Its geological unit is contributed to the west section of Tieling-Jingyu paleo-uplift of North China platform. The outcrops consist primarily of Mesozoic volcanic tuff, basalts, and weak interlayers (fractured zone). The underlying bedrock makes up of Archaean granitic gneiss [24]. The local topography has been excavated into step-like relief in the study area. The average slope angles of the open-pit mine reach to 28°~31° in the western part, 30.5° in the middle part, and 31.8° in the eastern part [22,25].
Around FWOPM, there exist several large enterprise groups, a large number of residential zones, and an industrial zone including many factories ( Figure 2). However, serious geological disasters, especially large-scale landslides, have become a direct threat to these residents and companies in Fushun city due to excessive and long-term coal mining. The first landslide with a historical record happened in 1927 [26]. According to incomplete statistics, more than 100 landslides have occurred up to now [25]. The main cause of landslides is primarily due to the stress imbalance resulted from high steep slope mining, which not only make the original slits cracked, but also a large number of new fissures appeared.
The native rock structures were changed, resulting in the decrement of rock strength, and making side slope rock deformation. Furthermore, the types of rock and lithology in the open-pit mine are highly variable, several large faults located under the open-pit, heavy rainfall and ground water infiltration would lead to landslides, too [27]. Therefore, it is significant to investigate the geological disaster situation to determine the potential landslides and to monitor large-scale slopes in motion before the slope loses steadiness and deep seated landslides occur.  The SAR data used in this study were acquired from the Phased Array type L-band Synthetic Aperture Radar (PALSAR) sensor onboard the Advanced Land Observing Satellite (ALOS). The PALSAR sensor operated in L-band frequency (1.27 GHz) with a ground surface resolution of 10 m × 10 m in HH polarization and a swath width of 70 km [28]. Table 1 lists the acquisition dates of the used SAR data. All scenes from the ascending orbit are in fine beam single polarization mode (FBS). Table 1. Summary of the used synthetic aperture radar (SAR) data provided by the phased array type L-band synthetic aperture radar (PALSAR) sensor onboard the advanced land observing satellite (ALOS).

Methodology
In order to obtain the displacement field time-series of large-scale landslides in FWOPM, a comprehensive data processing strategy was proposed, including the conventional DInSAR-SBAS data processing, MAI-SBAS data processing, and a combination of DInSAR-SBAS and MAI-SBAS results.
Conventional DInSAR-SBAS algorithm produces multiple unwrapped differential interferograms using SAR image pairs characterized by small baselines, which were used to generate the displacement time series via SVD (singular value decomposition) method. A detailed introduction of the SBAS algorithm can be seen in the study by Berardino et al. [29]. However, because of the low coherent vegetated regions, errors in the DEM, incomplete orbital errors removal, and residual atmospheric errors, it is possible that, at the end of the DInSAR-SBAS processing, some errors are still present. Therefore, we utilized two data processing steps to improve the accuracy of DInSAR-SBAS (bold contents in Figure 3). Firstly, a more reliable three-dimensional minimum cost flow (3D MCF) phase unwrapping algorithm consider the temporal information to help unwrap the low coherent areas by looking at the other coherent pairs applied to an irregular triangular network. Secondly, the L-band ALOS PALSAR data with significant residual topographic noise due to ALOS orbital parameters are not random in time [30]. Therefore, an additional method is used to remove the strong residual topographic noise.  The present method referred to as MAI-SBAS explored to obtain large azimuthal displacement time series by coupling the MAI and SBAS approach (Figure 4). First, we follow the same rule of the SBAS method to divide all the available SAR images into pairs characterized by small temporal and spatial baselines. Then, the selected image pairs were used to calculate the azimuthal displacement via a modified data processing of MAI to improve its accuracy. Finally, we apply the SBAS inversion strategy to retrieve the azimuth displacement time series.  The MAI method can obtain the azimuthal displacement through azimuth band splitting and interferometry technique [16]. According to MAI principles, the forward-and backward-looking images can be obtained from the raw SAR Single Look Complex (SLC) images by Azimuth Common Band Filtering (ACBF) method using modules of GAMMA software firstly [31]. For an interferometric pair, four new SLCs including the forward-and backward-looking master images, the forward-and backward-looking slave images can be produced. Then, the forward-and backward-looking pairs can be used to generate the forward-looking interferogram and the backward-looking interferogram. Subsequently, the MAI phase can be generated from the phase difference between forward-and backward-looking interferograms, which is proportional to the azimuth displacement. Therefore, the azimuth displacement dAZI can be calculated from MAI phase by the following equation [16].
where dAZI is the azimuth displacement, ΦMAI is the MAI phase, l represents the effective antenna length, n is the fraction of the full synthetic aperture width. Due to the azimuth band splitting, the MAI precision decreased significantly due to a higher phase noise and more processing artifacts than the conventional InSAR.
It is noteworthy that the forward-and backward-looking interferometric phases are differenced to produce the MAI phase. This step reduces some error sources including topographic and tropospheric contributions, which are common in conventional InSAR processing [16]. Nevertheless, the original MAI approach neglects the phase distortion from flat-earth and topographic effects due to the perpendicular baseline difference between forward-and backward-looking pairs. Jung et al. demonstrated that these residual flat-earth and topographic errors caused by the small difference between the perpendicular baselines of the forward-and backward-looking SLC pairs should be corrected before producing azimuth displacement. Theoretically, if the perpendicular baseline difference can be calculated accurately, the flat-earth and topographic phase residuals can be estimated and removed. However, the accurate estimation of perpendicular baseline difference is nearly impossible in the real MAI processing. Therefore, polynomial models were used to simulate and eliminate both phase residuals following the method described by Jung et al. [21]. For the displacement of the investigated landslide in FWOPM, the characteristic of deformation is that its motions occurred in both vertical and north-south directions. Therefore, the MAI results with north-south deformation are constrained by the vertical deformation and high coherence maps derived by D-InSAR method. In order to overcome the topographic, tropospheric and orbital noises during the data processing for L-band ALOS-PALSAR data, three steps were used to remove the topographic, tropospheric and orbital errors in this study. For both DInSAR-SBAS and MAI-SABS, the generation of interferometric phases can reduce some error sources including topographic, orbital and tropospheric contributions. Then, the polynomial model described by Jung et al. [21] was used to simulate and eliminate both residual flat-earth, orbital and topographic errors. Furthermore, the standard SBAS inversion steps include one step to perform orbital refinement and re-flattening process for all pairs.
In order to map the displacement field of slope, the combination of LOS displacement and azimuth displacement for ascending mode of observation is carried out. Let = [ ] be the three component vectors in north, east, and zenith, representing the displacement field, as shown in Figure 5, the ascending mode has the following relation [32]: sin sin sin cos cos

Data Processing and Results
The first step of data processing is the selection of image pairs with the proper perpendicular and temporal baseline that will be processed by DInSAR-SBAS and MAI-SBAS. The proper image pair selection is accomplished by imposing a threshold on the maxima perpendicular and temporal baseline, which means only small baseline interferograms will be generated. Considering the L-band ALOS PALSAR data with large perpendicular baselines are coherent, we set the perpendicular baseline to be 1500 m and the temporal baseline to be 3.5 years. The number of available interferograms is 22 and the more detailed information of the generated image pairs are listed in Table 2, including the acquisition dates of master and slave images, the perpendicular baseline, and the temporal and small baseline subsets. The distribution relevant to the image pairs and interferograms in the perpendicular and temporal baseline plane used in this study is shown in Figure 6.

LOS Displacement from DInSAR-SBAS
The standard two-pass differential InSAR approach was utilized to obtain the LOS displacements for all selected image pairs. As part of the DInSAR processing, an external digital elevation model (DEM) data from the Shuttle Radar Topography Mission-3 (SRTM-3) with 90 m horizontal resolution was used to remove the topographic component of interferometric phase. Before unwrapping the phase, it is necessary to improve fringe visibility and reduce phase noise. For this purpose, all the interferograms were processed using a multilook operation with six looks in the azimuth direction and four looks in the range direction to suppress the phase noises. Furthermore, an adaptive interferogram filtering algorithm was applied to the flattening image by setting large rectangular patch (window) sizes and high values of exponential power spectrum filter parameters [33]. In this study, we found that the 256 × 256 patch size and 1.0 filter parameter can improve the fringe visibility significantly. After that, the interferograms were unwrapped using the 3D minimum cost flow (MCF) algorithm dealing with Delaunay triangulation network [34,35]. Finally, the standard SBAS inversion steps are performed to obtain the LOS displacement of each time period. The LOS displacement can be converted into vertical and horizontal components. In this case, the GPS measurements indicated that the displacement in east-west direction is very small. In consideration of vertical component is a one dimensional vector, the conversion of the LOS displacement to vertical displacement can be expressed as [36]: where ddisp is the vertical displacement; dLOS is the LOS displacement; α is the incident angle, which can be obtained from the SAR data, the average incident angel is about 38.7° in this study. The geocoded vertical displacement mapped in WGS84 coordinate system is shown in Figure 7. It can be seen that landslides occurred at the edge of the huge open-pit with large area in the north part of FWOPM. The largest vertical displacement is estimated in Figure 7 to be up to 1.1 m.

Azimuthal Displacement from MAI-SBAS
The MAI approach was applied to obtain the azimuthal displacement for all selected image pairs in this section. Firstly, the forward-and backward-looking SLC scenes were produced using a normalized squint of n = 0.5. The same processing as DInSAR is then applied to generate the forward-and backward-looking interferograms, but the unwrapping operation was not required. Then, we conducted conjugate multiplication between the forward-and backward-looking interferograms to yield the final MAI interferogram. After that, the standard SBAS inversion steps are also performed to obtain the azimuth displacement of each time period. The geocoded north-south displacement maps in WGS84 coordinate system is shown in Figure 8. It can be seen that the region of horizontal deformation is similar to the area of vertical displacement, and the motion is direct to the bottom of the open-pit mine along the slope direction. The largest horizontal displacement is estimated in Figure 8 to be up to 1.45 m.

Figure 8.
North-south displacement maps in FWOPM obtained by MAI-SBAS approach. The red color represents the southward displacements for each time period, while the blue color represents the northward displacements for each time period. CP01 and CP02 represent two locations of GPS measurements. The background is a shaded relief map of SRTM DEM.

Integrated Displacement Velocity Field from DInSAR-SBAS and MAI-SBAS
In order to determine the deformation velocity field, the velocity results of DInSAR-SBAS and MAI-SBAS are integrated by Equation (2) described in Section 3. Figure 9 shows the integrated displacement velocity field map of FWOPM landslide for the period of January 2007-January 2011. The color represents the vertical displacement rate for the studied period. The red arrows represent the horizontal displacement velocity with different directions in the north and south parts of FWOPM.
It can be seen that the deformation velocity field of FWOPM landslide is measured with high spatial resolution. This resolution is fine enough to characterize the spatial variability pattern of the landslide body. The map indicates that the most active part of the landslide is the upper part in both north and south of FWOPM edges and that the deformation velocity decreases downslope. The maximum measured horizontal displacement rate is 487 mm/yr, while the maximum measured vertical velocity is 372 mm/yr.

Comparison with GPS Measurements
In order to quantify the DInSAR-SBAS and MAI-SBAS displacements derived in this study, we conducted a comparison of the InSAR displacement with the displacements measured on-site by two dual-frequency GPS receivers (the two locations CP01 and CP02 in Figures 7 and 8) at the north of FWOPM. The GPS data are processed by the GAMIT/GLOBK software [37]. Figure 10 presents the GPS displacement time series obtained for the periods January 2008-January 2011 for CP01 and November 2008-January 2011 for CP02 in black color. The corresponding amounts of the cumulated displacement measured from the InSAR method are depicted in red color. It can be seen that the deformation tendencies from InSAR were consistent with the deformation obtained by GPS. For the two components (vertical and north-south directions), the root mean square error (RMSE) between MAI-SBAS results and GPS measurements is in a range of ±10.5 cm for the north-south components, while the RMSE between the vertical displacements derived from DInSAR-SBAS results and GPS measurements is in a range of ±0.95 cm for the vertical component. Therefore, the landslide kinematic pattern is well identified by the proposed method.

Comparison with Optical Satellite Image
In this study, only two GPS stations located in the north part of FWOPM were available for comparison during the studied period. Therefore, it is impossible to perform a comprehensive verification between GPS measurements and InSAR measurements. In 2014, an accelerated landslide appeared in the south part of FWOPM. The clear ground fissures can be seen in Figure 11b in the south part of FWOPM caused by a landslide that is revealed by the high spatial resolution optical image from China's Gaofen-1 satellite in 2014, which provided a better chance to validate the InSAR results with the optical satellite image.
We conducted a comparison of the boundary of deformation region revealed by optical satellite image with high resolution and InSAR measurements. It can be seen that the detected deformation area by InSAR method in the south slope agree well with the location of the accelerated landslide region revealed by optical satellite image with high spatial resolution. In addition, the results in this study agree well with the results published by Nie et al. in 2014 [22]. It suggests that the deformation pattern can be found before the significant landslide occurred at least two years before. It also provides an important evidence for the landslide disaster controlling in the north part of the open-pit mine. Figure 11. The yellow dashed line marks the boundary of the south landslide in FWOPM [22]. (a) The color represents the integrated displacement rate for the studied period. The background is a shaded relief map of SRTM DEM; (b) The base image was captured by China's Gaofen-1 satellite in 2014, several ground fissures can be seen clearly in the south part of FWOPM; (c) The base image was an aerial picture of the study area from [22].

Conclusions
In this paper, we demonstrate the combination of DInSAR-SBAS and MAI-SBAS techniques applied to series of L-band SAR images for mapping and quantifying landslide displacement field. The typical north-south direction landslide at FWOPM was taken as an example. The 11 ALOS-PALSAR FBS data during the period from January 2007 to January 2011 were used to reach the aforementioned goals. The result in this paper reveals that the areas of south and north part of slopes have experienced significant northward/southward and vertical motions during the SAR data acquisition period, whereby the maximum accumulation in vertical direction exceeds 1100 mm and the maximum accumulation in north-south direction exceeds 1450 mm. The results have been validated by GPS measurements with a root mean square error of ±0.95 cm in vertical direction for DInSAR-SBAS and ±10.5 cm in north-south direction for MAI-SBAS, demonstrating that the proposed procedure was effective for monitoring large-scale horizontal and vertical displacements of landslides, especially for north-south direction landslides.
The proposed method is sensitive to the displacements along the vertical and north-south directions. However, it cannot provide comprehensive measurements for the motion along the east-west direction. Therefore, it is necessary to improve the accuracy of the offset-tracking based SBAS method to provide more displacement information. A coupling method among DInSAR-SBAS, MAI-SBAS, and offset-tracking-SBAS will be the subject of future investigations to monitor deformation in FWOPM and other places.