Analysis of Mining Waste Dump Site Stability Based on Multiple Remote Sensing Technologies

Abstract

The mining waste of open pit mines is usually piled-up in dump sites, making a man-made hill more than tens of meters high. Because of the loose structure of the dump sites, landslides or debris flow may occur after heavy rainfall, threatening local lives and properties. Therefore, dump stability analysis is crucial for ensuring local safety. In this paper, a collaborative stability analysis based on multiple remote sensing technologies was innovatively conducted at the Xudonggou dump of the Anqian iron mine. A small baseline subset (SBAS) analysis was used to derive the spatial and temporal distributions of displacements in the line-of sight (LOS) over the whole study area. The deformation in LOS is translated to the slope direction based on an assumption that displacements only occur parallel to the slope surface. Infrared Thermography (IRT) technology was used to detect weak aquifer layers located at the toe of possible landslide bodies. Then, numerical simulations based on the limit equilibrium method were conducted to calculate the factor of safety for three profiles located on the dump site. The results, emerging from multiple remote sensing technologies, were very consistent and, eventually, the landslide hazard zone of the Xudonggou dump site was outlined.

1. Introduction

Significant amounts of fragmented waste rock and loose top soil are excavated and piled at a certain site during the mining process, making up to several large man-made dumps with heights of more than one hundred meters [1]. The capacity of some waste dumps reaches tens of millions of cubic meters and some even reach hundreds of millions of cubic meters. Because of the loose structures of the dump sites, the unstable slopes pose a threat to local safety due to the high risk of landslide or debris flow, especially after heavy rainfall [2]. Therefore, stability monitoring of waste dumps is essential to ensure the safety of local lives and properties.

With the continuous advancement of science and technology, ground surveys such as prism leveling, total station measurements, and GPS measurements are widely used in China. However, it is difficult for ground survey technology to retrieve the deformation patterns of a whole area since it can only accurately measure deformation at a few points. As an important complement to traditional ground survey methods, synthetic aperture radar interferometry (InSAR) has become a popular tool for large area displacement monitoring, especially with the rapid development of time series InSAR (TS-InSAR) algorithms such as permanent scatterers interferometry (PSI), small baseline subset (SBAS) analysis [3,4,5,6,7,8,9,10,11,12,13,14,15,16]. Although, it has widespread applications and has been utilized to measure deformation with promising accuracy, it remains difficult to overcome the limitations of 1-D line-of-sight (LOS) measurements caused by the intrinsic side-looking geometry of SAR sensors [15]. Since ground displacements generally happen in both horizontal and vertical directions, how to translate the 1-D LOS measurements to the real displacements of the ground targets is a key problem plaguing the successful application of TS-InSAR technology.

At present, the ways of resolving three-dimensional surface displacements can be categorized in three ways [15,16,17,18,19,20,21,22,23,24,25,26,27]. The first way is by combining the multi-pass LOS and azimuth measurements generated from SAR observations. For example, three dimensional displacements can be extracted by combining multiple InSAR LOS measurements from at least three viewing angles [16,17,18,28,29], or by combining InSAR LOS measurements with azimuth measurements from offset tracking (OT) or multi-aperture interferometry (MAI). The second way is to integrate InSAR measurements with GPS data [30,31,32,33,34,35,36]. However, it is often difficult to satisfy the requirements of the above mentioned approaches. The above mentioned approaches all have their drawbacks and are only applicable in certain circumstances. For example, fused multi-pass InSAR LOS measurements can only be applied to areas where multiple stacks of data are available and high precision in the north-south directions is limited to high-latitude regions [18]. The fusion of InSAR LOS measurements with azimuth measurements obtained by OT or MAI has low accuracy in terms of north-south displacements because of the relative insensitivity to slow displacements of OT and MAI [17]. The integration of InSAR measurements with GPS data strongly depends on the number and distribution of GPS observations [35]. Fortunately, as a way to relax the requirements for TS-InSAR analyses, there is a third way that involves simply translating the InSAR measurements according to some prior information [37,38,39,40]. For example, the assumption that displacements only occur parallel to the ground surface can be made for landslides and glacial movements and horizontal deformation can be neglected when analyzing the subsidence of a large area, etc. Other remote sensing technologies, complementary to TS-InSAR, are also necessary for a complete investigation of ground displacement.

In this paper, a collaborative approach was proposed and conducted to evaluate the stability of the Xudonggou waste dump in the Anqian iron mine, using multiple remote sensing technologies. First, a TS-InSAR analysis is carried out to determine the LOS displacement of the Xudonggou dump site. Then, based on the assumption that ground displacement only happens parallel to the slope surface, the LOS measurements were translated to the slope direction. To assess the accuracy of the measurements translated from TS-InSAR, the retrieved displacements, along slope direction, were compared with the ground measurement data collected by total station instruments. In order to further investigate the deformation pattern of the dump site, aquifers on the slope were detected using thermographic imaging. Finally, the numerical simulation of dump stability was carried out based on two different assumptions. By using collaborative, multi-source, remote sensing technologies, a potential landslide danger zone could be precisely delineated.

2. Study Area and Dataset

2.1. Study Area

The Anqian iron mine is a large open pit mine, operated by the Ansteel Group Cooperation, with an annual total stripping capacity of 31.4 million tons, 13 million tons of iron ore and 8 million tons of ore processing. It is located in Anshan city, in the Liaoning province of China, which is marked with a red star in Figure 1a. There are three dump sites around this open pit, including Xudonggou, Yabaling and Xidabei. The mining waste is transported to the dump sites and piled up by trucks. Our study area, Xudonggou, located to the east of the open pit, is the oldest and most unstable dump site amongst them. As shown in Figure 1b, the dashed cyan line indicates the ridge line of the local topography and the mining wastes are generally piled up on the northern slope of the hill. Therefore, the structure of the Xudonggou waste dump can be categorized as having a side-hill fill structure. At dumps with this kind of structure, the instability is generally located on the north-facing slope, including the northeast and northwest. Therefore, the north slope of the Xudonggou dump site is generally an area of interest and where continuous monitoring is conducted.

Since the waste rock of the Anqian open pit mine has been dumped in Xudonggou for decades, the lack of stability of the north slope has become a threat to the safety of the mine as the dumping volume increases. It was closed due to the significant volume of dumping. The Xudonggou dump site reached a height of about 164 meters. In 2009, a landslide occurred at the dump. After the landslide event, two platforms were established along the slope, as depicted in Figure 1b. Platform 1 is on the lower part of the slope and has a width of approximately 25 meters. Platform 2 is on the upper part of the slope and has a width of about 10 meters. The dump’s slope gradient is approximately 32 degrees. Since two mining tailings transport pipelines and some local residents are located below the dump, landslides or debris flow in this area could be disastrous. Therefore, continuous monitoring of the Xudonggou waste dump is of great importance for the safety of local residents and their properties.

The dump site is mainly composed of a mixture of fragmented rock and loose top soil, see Figure 2. The rock comprises old metamorphic and mixed granites from the Archaea to the early Proterozoic period. Hard rock mass (e.g., low-grade magnetite quartzite, amphibolite, granulite and mixed granite) and soft rock mass (e.g., chlorite schist, quartz schist and phyllite) are both piled up in the dump. The hard rock bodies can withstand large shear forces without deformation, landslides rarely occur here. On the contrary, the soft rock bodies are generally plastic, which have low shear strength and deform easily.

In addition to the rock mass, fine-grained sand and silt are also found in this dump. According to their particle size, the sediments can be categorized into gravel (>2 mm) and sand (0.0625–2 mm). The sand can be further categorized into very coarse sand (1–2 mm), coarse sand (0.5–1 mm), medium sand (0.25–0.5 mm), fine sand (0.125–0.25 mm), and very fine sand (0.0625–0.125 mm). The materials with a particle size smaller than 0.0625 mm can be categorized as silt or mud [41]. Among these materials, the water permeability of silt and fine grained sand is relatively poor due to the small particle size. Therefore, water remains in these small particles after rainfall, resulting in aquifer layers in the dump site. Studies have shown that the lithology of soft rocks and aquifers is relatively weak; it is easy to form a weak soil- or mud-like layer under the influence of tectonic action, water, weathering and other external forces, and this could eventually result in a sliding surface or sliding zone [42,43].

2.2. Dataset

Previous studies have shown that most sliding body instability is caused by rainfall, in addition to individual rock and soil performance [44]. At the test site, the rainy season generally lasts from May until September. Therefore, the stability of the dump site during these four months is crucial to local safety. After looking at the archives, a descending stack of 13 COSMO-SkyMed images covering the study area and acquired from 6 May to 11 September 2015 was collected and employed in our time series analysis. As one of a new generation of high-resolution SAR satellite missions, COSMO-SkyMed is a military-civilian dual-use satellite constellation of four satellites established by the Italian Ministry of Defense and the Italian Space Agency [45]. The constellation works in an X-band, with multiple imaging modes and full polarization. The images are acquired by stripmap mode, with similar look angles of approximately 27 degrees. The thresholds of 0.27, 200 days and 680 m were applied for coherence, temporal and perpendicular baselines, respectively. As a result, 35 interferograms were generated from the 13 SAR images. The spatial and temporal baseline distributions are depicted in Figure 3, where the asterisks indicate the SAR images and blue lines represent the interferograms. The topographic phase was simulated and removed from the interferograms using the TanDEM-X DEM of 3-arc-second resolution (with spatial sampling of 90 m × 90 m) covering the study area [46,47]. TanDEM-X (TerraSAR-X add-on for Digital Elevation Measurements) is a German earth observation radar mission that consists of a SAR interferometer built by two almost identical satellites flying in close formation. With a typical spatial baselines between the satellites of 120-150m, a global Digital Elevation Model (DEM) has been generated. The absolute height error of TanDEM-X DEM is with about 1m an order of magnitude [47].

3. Methodology

3.1. Small Baseline Subset Analysis

Current TS-InSAR approaches are generally categorized into two groups. The first group uses interferograms generated in reference to a common master, e.g., the traditional PSI approach. The second group uses only good interferograms generated from all possible image pairs, e.g., the SBAS approach. Since the vegetation in mining areas has been seriously damaged, most of the area is composed of fragmented rock, gravel, sand and silt. The traditional PSI algorithm developed for urban areas does not work very well in the study area [5,6]. Besides, at least 20 SAR images are necessary in order to get a robust time series analysis result using traditional PSI. On the other hand, the small baseline subset method makes full use of the interferograms and, with small spatial and temporal baselines, can overcome the limitations of decorrelation with a reduced amount of SAR images [7,8,9]. Compared with the original SBAS method which only operates with isolated subsets of interferograms, the modified SBAS approach in the Stanford method for persistent scatterers (StaMPS) software can ensure temporal continuity by connecting separated subsets of interferograms with larger baselines [12,13]. It has been widely applied in the displacement monitoring of natural scenes with limited amounts of data [48]. Due to the above mentioned reasons, the modified SBAS approach in StaMPS developed for displacement monitoring of rural areas was well suited to this study [12].

The SBAS approach in StaMPS first generated interferograms with given thresholds on spatial and temporal baselines. It then connected separate subsets with larger baselines to ensure temporal continuity [12]. In the StaMPS SBAS, slow decorrelating filtered phase (SDFP) pixels were selected by carrying out amplitude dispersion and phase characteristics analysis. Then, the spatially correlated phase contribution was estimated and removed with bandpass filtering in the frequency domain. The spatially uncorrelated look angle error was then removed according to its correlation with the perpendicular baseline. After removing these two phase terms, a three dimensional phase unwrapping was carried out. The unwrapped phase of each pixel could then be divided into components based on surface displacement, topographic error, atmospheric disturbance, orbital error and noise. By applying iterative filtering in spatial-temporal domains, these different components could be successfully separated. Afterwards, time series deformation in the line-of sight (LOS) direction could be extracted using singular value decomposition (SVD) [12].

Due to the intrinsic, side-looking geometry of SAR satellites, the estimated displacement in LOS is just a projection of the real displacement. On the other hand, the deformation along the slope direction was measured by total station instruments. In order to assess the accuracy of the SBAS results, the deformation in LOS direction needed to be translated to the slope direction. For landslides, the simple assumption that the displacement was parallel to the slope surface was made [17,37,38,39]. Since the COSMO-SkyMed images were collected in descending orbit, and the sensor worked in a right side-looking direction. We will take Point A in Figure 4 as an example; the displacement along the slope direction can be split into vertical and horizontal components, as shown in the following equation:

$$\begin{array}{l}{d}_{A\_v}=d_A·\sin (\beta )\\ d_{A\_h}=d_A·\cos (\beta )\end{array}$$

(1)

where $d_A$ is the displacement of Point A parallel to the slope direction, β is the gradient of the slope, and $d_{A\_v}$ and $d_{A\_h}$ are the vertical and horizontal components of $d_A$ respectively. As shown in Figure 4, line OA represents the horizontal displacement $d_{A\_h}$ of Point A. Line CA represents the projection of $d_{A\_h}$ onto the vertical plane of LOS, which can be expressed with the following equation:

$$d_{A\_h-los}=d_{A\_h}·\sin (\pi -{\alpha }_A+{\alpha }_0)$$

(2)

where $d_{A\_h-los}$ is the projection of $d_{A\_h}$ onto the vertical plane of LOS, and ${\alpha }_A$ is the angle of line OA in reference to the south. ${\alpha }_0$ is the heading angle of the descending orbit with reference to south.

As depicted in Figure 5, the LOS displacement of Point A ($d_{A\_los}$) comprises the projections of $d_{A\_v}$ and $d_{A\_h-los}$ to the LOS direction, which can be expressed with the following equation:

$$d_{A\_los}=d_{A\_h-los}·\sin (\theta )+d_{A\_v}·\cos (\theta )$$

(3)

where $\theta$ is the look angle of the SAR satellite. By putting Equations (1) and (2) in Equation (3), the relationship between $d_A$ and $d_{A\_los}$ is expressed by the following equation:

$$d_{A\_los}=d_A·\cos (\beta )·\sin (\pi -{\alpha }_A+{\alpha }_0)·\sin (\theta )+d_A·\sin (\beta )·\cos (\theta )$$

(4)

Therefore, the deformation of Point A along the LOS direction is translated to deformation along the slope direction via the following equation:

$$d_A=\frac{d_{A\_los}}{\cos (\beta )·\sin (\pi -{\alpha }_A+{\alpha }_0)·\sin (\theta )+\sin (\beta )·\cos (\theta )}$$

(5)

Similarly, the displacement of Point B along the slope can be retrieved with the following equation:

$$d_B=\frac{d_{B\_los}}{\cos (\beta )·\sin (\pi -{\alpha }_B+{\alpha }_0)·\sin (\theta )+\sin (\beta )·\cos (\theta )}$$

(6)

where ${\alpha }_B$ is the angle of line OB with reference to the south direction. In the study area, ${\alpha }_0$ is approximately 11 degrees, $\theta$ is approximately 26 degrees, and β is approximately 32 degrees. The parameter ${\alpha }_A$ or ${\alpha }_B$ varies for different points due to their various locations on the dump slope.

3.2. Infrared Thermography

Infrared thermography (IRT) is a science dedicated to the acquisition and processing of thermal information from no-contact measurement devices [49]. According to the electromagnetic radiation theory, any object at a temperature above absolute zero (T > 0 K) emits infrared radiation [50]. Similar to a common camera which forms an image using visible light, a thermographic camera is able to collect infrared radiation and translate it into digital signals with very high temperature resolution in real time [51]. Then, the digital signal is processed based on the thermal contrast between various targets, which by their temperature and emissivity form a thermal image. IRT has been widely used in many applications, such as medicine, building inspection, and nondestructive testing [51,52,53,54]. As well as in the detection of heat distribution, IRT can be used in moisture detection [55].

Instead of physical temperature, the radiant temperature of an object was obtained by IRT technology. Generally, it is complicated to detect the physical temperature of an object using IRT. Fortunately, knowing the physical temperature of the object is not always necessary. It is sometimes enough to monitor the status of an object merely using the radiant temperature. A thermographic camera first collected the total radiation energy of the whole illuminated area. Then, the total radiation energy was used to calculate the radiant temperature of the object according to a function model. In fact, aside from collecting the radiation of the object, radiation from the background environment and atmosphere were also collected by the thermographic camera. Therefore, the radiation components from the surrounding environment and atmosphere had to be calibrated. The function model of the total radiation energy, collected by a thermographic camera with the radiant temperature of the object, is expressed by the following equation:

$$I_{mea}=I(T_{obj})\times \tau \times \epsilon +I(T_{sur})\times (1-\epsilon )\times \tau +I(T_{atm})\times (1-\tau )$$

(7)

where $I_{mea}$ is the total radiation energy collected by a thermographic camera, $I(T)$ is the radiation energy of a black body when the temperature is T, $T_{obj}$ is the temperature of the object, $T_{sur}$ is the surrounding temperature, $T_{atm}$ is the temperature of the atmosphere between the object and a thermographic camera, $\epsilon$ is the average emissivity in the spectrum of a thermographic camera, and $\tau$ is the average penetration rate of the thermographic camera in the spectrum [56]. Before collecting heat images, parameters such as emissivity, distance between the sensor and object, temperature of the atmosphere, relative humidity and temperature of the surrounding environment had to be inputted into the thermographic camera by the operator [57]. After that, the thermographic camera automatically translated the collected total radiation energy into the radiant temperature of the object according to the above mentioned function model.

3.3. Limit Equilibrium Method

Limit equilibrium (LE) methods have been widely utilized to obtain approximate solutions for stability problems in soil mechanics, due to their simplicity and accuracy [58]. There have been many LE methods developed for slope stability analysis, such as the ordinary method of slices, Bishop’s modified method, Janbu’s generalized procedure of slices and Spencer’s method [59,60,61,62,63]. These methods are based on the assumption that the failure surface is in simple shapes, e.g., a plane, circular arc or log spiral surface. These assumptions may not be particularly well founded, but often give acceptable results. Overall, Bishop’s simplified method is the most popular method since it is simple and yields results close to the rigorous limit equilibrium in many practical problems [64]. Therefore, the Bishop’s simplified method was used in the numerical simulation of the dump’s stability.

The slope of a dump site is more complicated than a natural slope in many ways. First, there are always several platforms along the slope, which makes mechanical analysis difficult. Second, the slope of a dump site is usually a mixture of loose top soil and fragmented rocks, instead of pure rock slope or soil slope. It is usually piled with gravel and the gaps are filled with soil brought by rainfall, which results in a special slope. Due to the complexity of the causal factors, the characteristics of the slope-forming materials are non-uniform and have a variety of severity, cohesion and internal friction angles. Because of the existence of multiple weak aquifer layers, the stability of weak layers in a dump site must also be attended to. In consideration of the variations in the parameters, stability analysis was carried out base on different parameters and a numerical simulation was carried out along three different slope profiles to estimate their LE state. The locations of the three profiles are depicted in Figure 6, where two platforms divide the MM′ and NN′ profiles into three parts and one platform divides the KK′ profile into two parts.

In the limit equilibrium analysis, the choice of parameters directly influences the results of the numerical simulation, no matter how the failure surface is assumed. The geometrical and mechanical characteristics of the slope were examined during field investigation. The values of these crucial parameters (e.g., the gravity density γ, cohesion and internal friction angle ϕ) were determined by statistical analysis on samples collected during the field investigation. It was generally believed that the internal friction angle of the slope depended on the mutual friction between different rock fragments because their particle grading and lithology were different. Since there is no cohesion between different rock fragments, the natural slope’s limit angle, formed when the mining wastes were first dumped, was considered the internal friction angle. After investigating the limit angles of 31 different spots on the dump site, an average internal friction angle (ϕ) of 38 degrees was selected for the simulation. The cohesion of the slope usually depended on the cohesive force of the soil between the gravel. The soil between the gravel is brought by rainfall, precipitated between the gaps and compacted after a very long time. Since the soil was already disturbed and reshaped, it was hard to collect the original samples to detect the cohesion. In addition, the soil was filled by rainfall and its characteristics changed greatly, which was identical to mixed surface filling soil. Therefore, a cohesion value of 12, for mixed filling soil in this area, was used for the simulation.

In the numerical simulation of the Xudonggou dump site, two types of failure surfaces were simulated, a plane surface and a circular arc surface. In the first simulation, a plane failure surface passing through the slope toe was assumed and the slope body above the failure surface was considered as a triangular prism, as shown in Figure 7. In consideration of cohesion, the factor of safty (FOS) was calculated according to the following equation,

$$K_2=\frac{tg(\varphi )}{tg(\delta )}+\frac{2·c·\sin (\beta )}{\gamma ·H·\sin (\beta -\delta )·\sin (\delta )}$$

(8)

where ϕ is the internal friction angle, β is the gradient of the slope, H is the height of the dump, and δ is the slope angle of the failure surface. In this simulation, a dangerous failure was determined by iteratively finding a certain slope angle δ between 0 and β, which made K₂ the smallest. The FOS of each part of the profile was calculated first, then a total FOS was calculated for each slope.

Taking the cohesive force of soil into consideration, the failure surface was assumed in the second simulation to be a circular arc. By dividing the landslide body into small soil strips, the factor of safety could be calculated with the following equation:

$$F_s=\frac{{\displaystyle \sum K_i}}{{\displaystyle \sum (W_i+Q_i)·\sin {\alpha }_i}}$$

(9)

where $W_i$ is the gravitational force of the ith soil strip, ${\alpha }_i$ is the slip angle at the bottom of the ith soil strip, $Q_i$ is the vertical force on the i th soil strip. $K_i$ and $m_{\alpha i}$ can be calculated with the following equations:

$$K_i=\frac{c_{ti}{b}_i+U(W_{ti}+Q_i)\tan {\phi }_{ti}}{m_{\alpha i}}$$

(10)

$$m_{\alpha i}=\cos {\alpha }_i+\frac{\sin {\alpha }_i\tan {\phi }_i}{F_s}$$

(11)

The meanings of the symbols in Equations (9–11) refer to the specifications for design of highway subgrades (JTG D30-2015) of China. The FOS in Equation (9) was obtained by summing all the soil strips respectively. In this simulation, a comprehensive search was carried out according to different arc centers and radius in order to find the smallest FOS of each profile. The FOS of all the profiles was calculated respectively, as well as the overall stability factor.

4. Discussion of Experimental Results

4.1. SBAS Results

The displacement velocity map along the line of sight (LOS) direction over the Xudonggou dump site is depicted in Figure 8. The Google Earth image was used as background image. As shown in Figure 8, negative displacement rates indicated projected movements away from the satellite platform, whereas positive values indicated projected movements toward the satellite platform. In the study area, a total of 3072 pixels were selected for a time series analysis. Most of the points were located at the platforms or roads and the highest density of PS points was on the south slope. However, only two groups of PS points were selected for the area in which we were interested. These are highlighted by the red and yellow rectangles.

As shown in Figure 9, The LOS displacement velocity of PS points basically obeyed the normal distribution. The mean displacement rate over the study area was about −1.44 mm/yr with standard deviation of 8.22 mm/yr and approximately 73% of the pixels showed displacement velocities between −9.66 (−1.44 − 8.22) mm/yr and 6.78 (−1.44 + 8.22) mm/yr. In statistics, the standard deviation is a measurement used to quantify the dispersion of a set of data values. In the SBAS results, a lower standard deviation indicated that the displacement rates of PS points tended to be close to the mean displacement rate, whereas a higher standard deviation indicated that the displacement rates of the PS points were spread out over a wider range. Generally speaking, a lower standard deviation value is expected in time series InSAR results. This is because most PS points should be more or less stable in a study area and deformation only happens to a small portion of the PS points. The largest displacement velocity along the LOS was −42.17 mm/yr, whereas the lowest was about 0 mm/yr.

For landslide studies, even if the assumption that deformation only occurs parallel to the slope surface is applied for simplicity, the true deformation direction of PS points in rugged areas is still a complicated issue due to the diversity of local incidence angles. The deformation estimated by SBAS was just a projection of the true deformation onto the LOS direction. If a very high precision DEM is available, the three dimensional projection geometries for all PS points can be automatically derived according to their respective geographic locations. Then, accurate translation from the LOS direction to the slope direction can easily be carried out. Unfortunately, there is no accurate DEM for the Xudonggou dump site to carry out automatic translation for the 3072 PS points. Therefore, geometric translation was only conducted for four points near the terrestrial measurement points for the purpose of cross validation.

From May to September 2015, total station measurements were conducted weekly at the study site. The locations of total station measurement points (A1-A8 and B1-B5) are depicted in Figure 1b. The cumulative deformation of the ground measurement points along the slope direction are shown in Figure 10, indicating that the cumulative displacements of the 13 points varied from 29.06 mm to 137.94 mm. Along Platform 1, A8 presented the largest movement (about 117.73 mm), whereas A3 presented the smallest deformation (about 46.47 mm). Along Platform 2, B1, B4 and B5 were relatively stable with displacements of 32.87 mm, 51.72 mm and 29.06 mm, respectively. On the other hand, B2 and B3 showed relatively poor stability, with cumulative deformations of about 130.81 mm and 137.94 mm respectively.

The deformation of the nearest PS point along the slope was compared with the ground measurement data of Points A3, B1, B2, and B3, respectively. For convenience, the cumulative deformation from SBAS was resampled according to the time lines of the ground measurements. As shown in Figure 11, the cumulative displacements from SBAS and ground measurements have been consistent over all the points, especially from May to the end of July.

In August and September, SBAS underestimated the displacements of B2 and B3, highlighted with a green box in Figure 11c,d. This is probably because SBAS was only able to measure displacements correctly within the wrapped phase (−π,π), which was approximately 31.2 mm in the LOS direction according to the wavelength of the X-band SAR sensors. By translating the measurement extent of the LOS to the slope direction, according to the specific locations of Points B2 and B3, a measurement extent of 43.33 mm was obtained. However, the cumulative deformation of Points B2 and B3 were far beyond the measurement limits of the X-band SAR, which are 130.81 m and 137.94 mm, respectively. Therefore, the cumulative displacements for Points B2 and B3 were underestimated in August and September. On the other hand, fortunately, there was no underestimation for Point A3 due to its smaller cumulative displacement.

The authors also plotted the weather from 6 May to 15 September in Figure 12, where the yellow bars represent a sunny day, the green bars represent cloud, the cyan bars represent light rain or showers, and the red bars represent moderate rain or heavy rain. As is clearly shown in Figure 12, moderate continuous and heavy rainfall happened in August, which explains the dramatic increase in displacements in the ground measurements. The continuous rain may be another reason for the underestimates of SBAS in August, because of the short wavelengths of the X-band SAR sensors.

For the purpose of quantitatively assessing the performance of SBAS, the relative error and root mean square error (RMSE) of SBAS measurements, with reference to the ground measurements of the same time, were calculated. Let us take Point A in Figure 4 as an example, the relative error, $erro{r}^{\prime }$, and RMSE of Point A can be expressed with the following equations:

$$\begin{array}{l}erro{r}^{\prime }=\frac{{\displaystyle \sum \left(d_{A\_SBAS}\left(t_i\right)-d_{A\_Total\_Stations}\left(t_i\right)\right)}}{N}\\ RMSE=\sqrt{\frac{{\displaystyle {\sum }_{i=1}^N{\left(d_{A\_SBAS}\left(t_i\right)-d_{A\_Total\_Stations}\left(t_i\right)\right)}^2}}{N}}\end{array}$$

(12)

where $d_{A\_SBAS}\left(t_i\right)$ is the displacement along slope direction estimated by SBAS at moment $t_i$, $d_{A\_Total\_Stations}\left(t_i\right)$ represents the displacement of ground measurements at moment $t_i$, and N is the number of observations conducted during the temporal span. RMSE is a measure of accuracy, which is the square root of the average squared relative errors. In general, a lower RMSE indicates that the SBAS result is more consistent with the terrestrial measurements. The estimated relative error and RMSE of the four points over the whole period is shown in

Table 1. Relative Error and Standard Deviations.

Point ID	A3	B1	B2/B2_Corrected	B3/B3_Corrected
$erro{r}^{\prime }$ (mm)	−2.26	8.94	−17.18/−7.89	−20.99/−6.20
RMSE (mm)	8.09	4.79	14.54/5.91	22.66/12.46

. Generally speaking, the accuracy of the SBAS result for Point A3 was the highest according to the relative error index, whereas the worst was for Point B3. If all the measurements from May to September are taken into consideration, the relative errors and standard deviations of Points B2 and B3 are even higher than that of Point B1. However, if the obviously erroneous estimates in August and September are neglected, the relative errors and standard deviations of Points B2 and B3 are much smaller. By comparing the RMSE index, it is obvious that Point B1 has the smallest residual error, whereas Point B3 has the largest. If the clearly biased estimations in August and September are discarded, the RMSE of Point B2 can be significantly improved.

4.2. IRT Results

In this study, IRT technology was used to detect aquifers in the dump site. Since the granularity of sand is smaller than that of rocks, it usually forms water-bearing aquifers. Due to the large thermal capacity and thermal inertia of water, the temperature change of the aquifer was less than that of dry rock. The evaporation of water also made the temperature of the aquifer decrease. Therefore, the radiant temperature of the aquifer was lower than that of the surrounding rocks in the thermal infrared images. According to the above principle, aquifer layers formed by water-containing sand could be detected using IRT.

The thermographic camera used in this study was a VarioSCAN-3021ST. It has a very high thermal resolution of 0.03 K with a working wavelength of 8~12 um. The calibration range varies from −40 °C to 1200 °C. In addition, a high-resolution digital camera (Cannon 40D) was used to shoot the study area. By comparing the thermal images and the photos, aquifers could be precisely located. The processed thermal image is shown in Figure 13, where we can see the temperature of the three areas (highlighted in blue rectangles) was lower than that of the surrounding areas, with a difference of approximately 5–6 degrees. All three low-temperature areas were sand layers according to the field investigation, indicating aquifers.

The instability of a sliding body is closely related to inherent mechanical properties and external triggering factors [58]. The mechanical properties included the slope gradient and the characteristics of the slope-forming materials. The external triggering factors included rainfall, snow melt, ground water infiltration, volcanic eruption, seismic events and human activities. Previous studies have shown that most sliding body instability is caused by rainfall, in addition to its own rock and soil performance [44]. Rainfall infiltration increases the pore water pressure on the sliding material, thereby reducing its own shear strength and causing instability of the sliding body. The shear strength of the sand layers is relatively weak. Due to the relatively small granular size of sand layers, the pore water pressure even increases its ability to slide after intense rainfall [59]. According to engineering geology, aquifers are usually found at the lower part of a large sliding body [44]. Therefore, the slope body above the aquifer layers is the dangerous region which needs to be continuously monitored.

4.3. Numerical Simulation Based on the LE Method

In LE simulation, the factor of safety (FOS), which is the ratio of anti-sliding force to sliding force, is calculated for each slope respectively. Theoretically, a FOS less than 1.0 indicates that there is already sliding damage. The specifications for design of highway subgrades (JTG D30-2015) stipulates that the FOS threshold for a road slope is 1.35, meaning that coefficients smaller than 1.35 indicate instability. Unfortunately, there is no clear regulation for mining dump sites. Taking some unforeseen complex situations into consideration (e.g., abnormal weather phenomena), a FOS larger than 1.2 is considered to meet the requirements for safety and stability.

The simulation results for both assumptions are shown in Figure 14, where (a–c) are simulation results based on the plane sliding surface assumption, and (d–f) are simulation results based on the circular arc sliding surface assumption. The blue lines in Figure 14 represent the slope surface, and the red lines indicate the detected dangerous sliding surfaces for each profile. The simulated FOS in Simulation 1 is shown in

Table 2. FOS of each profile when the failure surface is a plane.

Profile Name	Slope 1 (Lower Part)		Slope 2 (Middle Part)		Slope 3 (Higher Part)		Total
	K2	δ	K2	Δ	K2	δ	K2	δ
KK′	2.70	21	1.69	29	-	-	1.69	29
NN′	2.71	19	1.51	36	1.56	32	1.08	36.5
MM′	2.46	21	1.74	28	1.68	29	1.49	28

. It is clearly shown that the total FOS of profile KK′ was the largest, whereas that of NN′ was the smallest of the three. A total FOS of 1.08 indicated that profile NN′ was the most dangerous slope in this simulation. By comparing the different coefficients of the three sections of profile MM′, we could tell that the higher section was more dangerous than the middle section, and the lower section was the most stable. The situation at profile NN′ was slightly different than that of MM′; the middle section was the most dangerous. In this simulation, only profile NN′ showed a FOS value less than 1.2. The other two profiles were generally stable according to Simulation 1.

The least favorable situation and potential failure surface was found in the second simulation, as shown in

Table 3. FOS of each profile when the failure surface is circular arc.

Profile	KK′	NN′	MM′
Most dangerous	1.4115	1.0139	1.1673
No. 2 dangerous	1.4921	1.2104	1.3492
Weak layer	-	1.4078	1.5643

. Compared to the first simulation, the FOS in the second simulation was generally smaller. Nevertheless, the FOS for profile NN′ was still the smallest, whereas the coefficient for profile KK′ was the largest of the three. The FOS of profiles MM′ and NN′ were both smaller than 1.2, indicating that the two profiles are quite unstable. The critical FOS value of 1.0139 of profile NN′ means that this profile has most likely moved.

The results from both simulations are consistent with each other. From the FOS values we can tell that profile KK′ is quite stable, with two FOS larger than 1.2. The smallest FOS of profile NN′ is 1.0139, indicating that this slope profile has probably been destructed. The obvious collapse and fissure found near profile NN′ during field investigation (Figure 15) also supports the conclusions of both simulations. Unfortunately, the FOS for profile MM′ is also smaller than 1.2, indicating that it is also quite dangerous.

5. Conclusions

In this paper, a collaborative analysis of the stability of the Xudonggou dump innovatively combined SBAS, IRT, total station measurements and numerical simulation based on the LE method. These technologies all have their advantages and drawbacks. SBAS is able to detect the spatial distribution and temporal variation of displacements over the whole study area with acceptable precision. Based on a simple assumption that landslides only occur parallel to the slope direction, the LOS measurements could be easily translated to the slope direction. Unfortunately, displacements orthogonal to the LOS direction could not be detected by SBAS due to the intrinsic side-looking geometry of the SAR satellites. Additionally, the phase unwrapping may continue to bias the results if the cumulative displacement exceeds the 2π measurement extent of the X-band SAR. Ground measurements using traditional instruments like total stations and prism leveling can provide displacements with very high accuracy. However, it is time consuming and requires significant manpower and resources. By merely collecting, measuring and comparing a limited number of points to SBAS, the higher cost of manpower and resources can be avoided. As a widely used technology in heat and moisture information collection, IRT is able to find water-bearing aquifer layers. These aquifers are probably weak layers due to their relatively small granular size, especially after intense rainfall. According to engineering geology, the slope body above aquifer layers is usually considered a dangerous area in terms of landslides. Numerical simulation based on limit equilibrium is a classic tool of engineering geology used to carry out slope stability analyses. It is easy to implement and often gives acceptable results to many practical problems, although sometimes the assumptions may not be well founded. By collaboratively using multiple remote sensing technologies, a full understanding of the Xudonggou dump site has been achieved. Based on the various measurements and numerical simulations, we are confident that the danger zones of the Xudonggou dump site have been recognized. They are shown in Figure 16. The collaborative use of multiple remote sensing technologies, as in this article, provides new insight into the various deformations caused by both human activities and natural disasters.