Carajás Mineral Province, located in the Brazilian Amazon region, encompasses the world’s largest iron reserves with mining exploration carried out through open pit benching. Open pit operations usually have significant areas of extent with rock mass movements and surface displacements that potentially lead to slope instabilities with risks to personnel and equipment and affecting production. It is worth noting that small surface movements on a mine highwall may be a signal of a failure. Thus, a key concern for the mining industry is the prediction of mining-induced deformations of ground surface.
Taking advantage of multitemporal SAR acquisitions, the use of advanced differential SAR interferometry (A-DinSAR) techniques improves the capability for detecting deformation phenomena. A differential SAR interferometry (DInSAR) time-series (DTS) technique has been proposed and successfully used22.214.171.124.–6 and consolidated using the small baseline subset (SBAS) concept.78.9.10.–11 A persistent scatterer interferometry (PSI) approach1213.14.–15 based on a stack of differential interferograms relies on identifying pixels whose scattering properties vary little with time and look angle, allowing a temporal analysis of the interferometric phase of individual point targets and providing accurate information related to the surface target displacements. PSI provides better spatial resolution and accuracy and can better model and diminish the effect of the atmospheric phase, which is related to the path delay heterogeneity during the acquisitions times, whereas the DTS technique can detect larger deformation rates and provides more homogeneous and extensive information related to the ground deformation, but at the cost of a loss in spatial resolution.
This paper is an outgrowth of the previous DInSAR investigations carried out with open pit iron mines in the Brazilian Amazon region.1617.18.19.–20 This study presents an application of A-DInSAR techniques for monitoring surface deformation in a large area of extent with several open pit iron mines located in Carajás Mineral Province (Brazilian Amazon region). In the next section, a description of this particular mining area, as well as some related ground deformation instabilities, is presented. The combined approach based on DTS and PSI techniques using 33 StripMap TSX-1 scenes is presented in Sec. 3. The results of the combined techniques are presented in Sec. 4. The potential and the results of this approach for monitoring ground movement in a large mining area are discussed in Sec. 5. In Sec. 6, the validation is presented with topographic measurements (total station/prisms). The conclusion is given in Sec. 6.
The most important Brazilian mineral province is located in the Amazon Region, Carajás-PA, encompassing 39 iron bodies with reserves of 18 billion tons in an area of . This province is marked by mountainous terrain, characterized by a set of hills and plateaus (altitudes from 500 to 900 m) surrounded by southern and northern lowlands (altitudes around 200 m), deep chemical weathering that produces thick oxisols (latosols), totally covered by ombrophilous equatorial forest communities with complex and multilevel canopies and numerous species.21
The iron deposits are covered by thick, hard iron-crust (lateritic duricrusts) developed over volcanic rocks and ironstones. Specific low-density savanna-type vegetation (campus rupestres) is associated with the deposits and shows a strong contrast (clearing) with the dense equatorial forest. Fully owned by Vale S.A. mining company, the exploration in Carajás is carried out through state-of-the-art open pit benching. The current mining activities are related to two iron ore bodies (N4 and N5) and started with the N4E mine in 1984 and with the N5W mine in 1998. Carajás Mining Complex (N4 and N5 bodies) (Fig. 1) produces metric tons of iron ore per day.17
Surface instabilities can be expected at any mining activity. Open pit operations usually have significant areas of extent and can also influence large portions of terrain adjacent to the pit crest. Rock mass movements and surface deformations potentially lead to slope instabilities or wall failures due to regular open pit mining operations. This scenario in Carajás gets worse over time due to intense deep excavations in saprolitic soil and rock masses of low geomechanical quality, coupled with blasting practices and heavy precipitation of the moist tropics, which has deleterious effects on the overall stability. Vale’s geotechnical team has monitored the presence of fractures on bench walls and tension cracks on berms and road ramps through visual inspection, total station/reflecting prisms measurements, and ground-based radar.
A stack with 33 TSX-1 StripMap images (repeat cycle of 11 days) was used in this investigation. Conflicts in acquisition programming of the satellite caused four interruptions in the original 1-year acquisition coverage starting on March 20, 2012, with four interruptions in 2012 (December 31) and 2013 (February 13, February 24, and March 7). The single look complex (SLC) images were acquired under ascending passes (), with an incidence angle range of 39.89 to 42.21 deg, a spatial resolution of (), pixel spacing of (), and a width swath of 30 km. To minimize the topography phase error in the interferometric process, a high-resolution DEM was produced based on a panchromatic GeoEye-1 stereo pair. The GeoEye stereo images were acquired over the study area on July 1, 2012. The first scene was collected with nominal azimuth and elevation angles of 29.4 deg and 82.4 deg, respectively, whereas the second scene was acquired with azimuth and elevation angles of 187.42 deg and 62.20 deg, respectively. The images were provided with 0.5-m spatial resolution and with rational polynomial coefficients. The generation of the DEM was based on OrthoEngine PCI Geomatics through the rational function method as the geometric model. The panchromatic DEM was produced using 2 m spacing, and its elevation values compared to seven well-defined accurate vertical check points have provided RMS and maximum errors of 1.2 and 1.6 m, respectively.22
The proposed approach described in this paper is shown in Fig. 2, where the combination of two techniques, DTS and PSI, are used to improve the capability of the PSI technique for detecting nonlinear deformation, removal of the orbit phase trend error, and removal of the atmospheric phase delay in a large area, as well as increasing the numbers of point density of the final results. The deformation model and the DEM errors estimated by the DTS are subtracted from the master referenced differential interferogram, and the residual phases are processed with PSI technique. This approach was already applied to detecting nonlinear deformation in a small area of the Carajás Mineral Province.19 In this work, we focus on the removal of the atmospheric phase component in a large area, where it is expected that it is not constant over the area.
The use of the standard DInSAR technique for monitoring surface deformation has been applied since the early ’90s. DInSAR technique aims to measure deformation on the ground using a pair of SAR images acquired at different times and different positions of the satellite between two acquisitions. The interferogram generated from these two coregistered images has phase’s component contributions from the topography, deformations, atmosphere, and noise. By knowing the positions of the satellite and the topographic surface, it is possible to subtract the topographic phase component and to measure the deformation. After removal of the phase component related to topography what remains in the differential interferogram is a contribution of the ground displacement between acquisitions added with other undesirable component, represented by
The ambiguity of in the differential interferometric phase (fringes) makes it impossible to interpret this phase in terms of absolute range change. However, it is possible to estimate the relative range change between points within a differential interferogram by integrating the number of fringes between them, using a process known as phase unwrapping.23 Ground deformation monitoring is suitable when the exceeds the other four phase components represented in Eq. (1). For this reason, the use of a precise DEM is desirable to compensate for part of the topographic phase error. A spatial filtering (multilook) can reduce the phase noise; however, the atmospheric phase delay and the residual phase due to orbit errors cannot be filtered out or estimated, respectively.
Taking advantage of multitemporal SAR acquisitions, the use of A-DinSAR techniques has improved the capability for detecting deformation phenomena. DTS and PSI algorithms have been developed to better address the standard DInSAR limitations using a redundant number of differential interferograms with the potential to determine spatially and temporally ground displacement, where the desirable deformation information can be separated from topography error, atmospheric delay, and noise.
Consider a set of differential interferograms generated from a set of SAR images acquired at the ordered time , where each differential interferometric pair is constructed in a given time interval (), as represented in Eq. (1), following the rule of small time interval between acquisitions or SBAS.7 The observed multilook unwrapped phase values of this point in relation to the reference point may be organized in a vector of elements, represented by Eq. (2).2 A physically sound solution can be found in terms of mean phase velocity among time-adjacent acquisitions using the SVD decomposition,7 and a final integration can be used to achieve the final solution .
The estimation of the topographic error () can be determined by the coefficient of a linear regression that relates the topographic phase with the perpendicular baselines that compose the interferograms, based on the following relation:
In this study, the inversion of the system of equation was performed using an extension of the SVD with a set of additional weighted constrain on the acceleration of the displacement to control the smoothness of the time-series solutions.3 Temporal smoothing is enhanced using a finite difference approximation constraint, assuming that less deformation occurs during a short period of time. The smoothing constraint and the height error-related term were incorporated into the inversion of the system that relates the observed unwrapped interferograms phase with the average displacement velocity.24
A stack with 33 TSX-1 StripMap scenes was used to perform this investigation. It was planned to have one image every 11 days, but this sequence failed twice, creating discontinuities of 22 and 44 days. To cover the time span of acquisitions, it simulated the differential interferometric pairs for a time interval up to 45 days and a maximum perpendicular baseline up to 400 m (Fig. 3). Visual inspection was used to discard interferograms with compromising phase unwrapping errors. Most of the selected pairs belong to the shortest time interval among acquisitions (11 days).
The results of the DTS, ground displacement () and topographic error () are resampled for full resolution to have the same pixel spacing as that of the differential interferograms for the PSI analysis (Fig. 2). An important transformation is performed to change the phase displacement () to the PSI phase model, which is referenced to the master image, as used in the PSI analysis, represented by
The PSI technique relies on identifying pixels whose scattering properties vary little with time and look angle in a stack of coregistered SLC images. The idea of PSI is to analyze the temporal and spatial characteristics of the interferometric phase of individual point targets. These point targets are coherent even for the interferometric pairs with long spatial baselines and remain stable over long time periods to permit analysis of the phase history.1213.14.–15
PSI analysis was carried out using the interferometric point target analysis (IPTA)14 software, which is the implementation of PSI by GAMMA remote sensing and consulting AG (GAMMA). The IPTA software is a toolbox that can support many different methodologies including different alternatives for PS candidate selection, spatial and temporal phase unwrapping, and supporting approaches for single as well as multireferenced stacks; see Eq. (11). The processing sequence included SAR SLC image coregistration to generate the stack of interferograms (the master scene was selected based on a configuration that has low perpendicular baseline dispersion and is near to the center of the image sequence), the stack of differential interferograms generation, and the point target candidate determination (based on spectral phase diversity and low intensity variability). Figure 4 shows the flow diagram outlining the processing sequence for the PSI analysis using the IPTA technique shown in Fig. 3. The IPTA is based on a two-dimensional (2-D) regression (phase model indicating a linear dependence of the phase displacement with time and a linear dependence of the topographic phase with the perpendicular baseline). For the first regression to find the displacement rate, the best fit is found within the predefined parameters (maximum and minimum rate intervals) and the standard deviation bound of the regression; for the second regression, related to the estimation of the topographic error, the best fit is searched within the predefined height error interval. The regressions provide a set of layer (stack) of the unwrapped phase () and residual phase (), as well as single layers of the topographic phase error (), the displacement rate (), and the standard deviation of the linear regression fit of the displacement (). The model refinement shown in Fig. 4 presents several processing steps: phase unwrapping correction of layers with patch effect, baseline refinement, atmospheric phase removal from the spatially filtered residue, updating of the PS interferometric phase from estimated displacement rate, and the topographic phase error. In the last step, the residual atmospheric phase layers are checked for a new iteration. The computation of the linear and nonlinear displacement and the corrected digital elevation model are carried out after the last iteration.
In this analysis, 33 TSX-1 images acquired from March 2012 to April 2013 were used. The configuration of the interferometric pairs in relation to the master image (October 4, 2012) is shown in Fig. 5. PS candidates were estimated based on the amplitude dispersion index and low spectral diversity at each pixel of the stack of the coregistered images.13,14 Considering the set of SAR images acquired at the ordered time , a reference point located in a stable area, and a selected generic PS, the observed wrapped phase values of this point in relation to the reference point can be organized in a vector of elements as
The results obtained from the DTS analysis, the topographic phase error , and the phase deformation model , represented by Eq. (6), are subtracted from the observed phase Eq. (7) to reduce the amount of these phase components, aiming at improving the PSI analysis with IPTA. The phase is subtracted (module-) from the observed phase vector Eq. (7), resulting in the first residual wrapped phase, which is then subtracted (module-) from the phase displacement model, as shown in Fig. 2, resulting in a residual wrapped phase represented in Eq. (8), which is used as input phase for the PSI analysis represented in Fig. 4.
Performing PSI analysis for a large area is always a challenging task regarding the removal of the atmospheric phase component. An important aspect of the IPTA technique is the possibility of a stepwise, iterative improvement for different parameters to estimate the topographic error, atmospheric phase delay, and ground displacement. The residue layers of the linear regression related to the displacement can present phase unwrapping errors, including the patch effect due to the local reference used. To overcome this problem, the layers with phase unwrapping error are transformed to a complex format filtered spatially and afterward unwrapped again until all layers have been unwrapped properly. This processing step can be repeated until a suitable box size for the filter is found.
In this work, three iterations were performed to improve the estimation of the ground displacement and the removal of the atmospheric phase components. A new iteration is decided based on the statistics of the residual atmospheric phase of all layers. The atmospheric phase component may account for the most of the linear regression deviation (residues) related to the linear deformation, considering that the larger phase components related to the linear and nonlinear movement and topographic error have been removed previously from the DTS results. The spatial statistics of the atmospheric residual phase for all layers presented a range of values from 0.5 up to 5.1 rad, 0.3 up to 0.8 rad, and 0.3 up to 0.7 rad for the three iterations, respectively. Figure 6 shows the estimated atmospheric phase component (in module) for one layer composed by the images master (October 4, 2012) and slave (October 28, 2012) in three iterations. It can be noted in Fig. 6 that, after the second iteration, the residual phase remains very low.
The final solution for the deformation was obtained by adding the PSI result to the DTS phase displacement model, represented by Eq. (9), and for the final topographic error by the addition of the components from PSI analysis with the DTS, as represented in Eq. (10). Only the PS present in the PS_mask (Fig. 2) was validated.
Results and Discussion
The ground displacement map obtained with the proposed methodology, using 33 TSX-1 scenes and covering the time span from March 2012 to April 2013 (dry and wet seasons), is shown in Fig. 7. PSI processing, using previous information about the phase displacement model and DEM error derived from the DTS, allows the detection of ground movement in a large area, covering the complex of the N4 and N5 open pit iron mines (Fig. 1).
No relevant ground displacement was found around and within the pits and mining infrastructures (green-bluish regions, Fig. 7). However, high deformation rates (yellow-reddish regions, Fig. 7) with negative values corresponding to motion away from the satellite were detected over the waste piles corresponding to the letters A, B and C (Fig. 7). Furthermore, it was also possible to detect evidence of deformation over the slopes of the N5W mine (letter D, Fig. 7). For the waste piles, the detected displacement was interpreted as related to settlements, showing values normally expected for this manmade structure. On the other hand, lithostructural and lithogeomechanical attributes played a key role in the displacements related to the cut slopes in N5W. The maximum accumulated displacements (along the LoS and away from the satellite) were (waste pile NW-1, letter A), (waste pile W, letter B), (waste pile S-IV, letter C), and (cut slopes of N5W, letter D).
Differential interferometric measurements shown in Fig. 7 were carried out with respect to a reference point that was selected within the area of interest, assumed as stable, and time referenced to the acquisition date of the first scene. Therefore, displacements data provided for PS are relative, not absolute, data. The accuracy of these measurements can be estimated through the dispersion of PS values regarding the reference point, expressed by the standard deviation of the displacement rate, represented according to Gamma25 by
Figure 8 shows the standard deviation map of the average displacement rates for the 33 images processed using the proposed methodology.
Geomechanical classes and the corresponding deformation rates on the ore bodies (Fig. 1) are shown in Fig. 9 for the N4W and N4E mines. It can be noticed from Fig. 9, for both mines, that the more relevant ground displacements are related to the poor and very poor geomechanical classes. This shows that, apart from the mining activity, poor geomechanical characteristics have played a key role in the displacement process. Table 1 reveals some relationships of PS and geomechanical classes for the mines N4W and N4E, showing that the highest amount of detected PS was associated with geomechanical classes with the largest areas but presenting the lowest densities (classes IV and V). These two classes present geomechanical attributes of poor quality and the most intense exploitation activities.
Relationship between PS and geomechanical classes for the N4W and N4E mines.
|Geomechanical class||Area (km2)||Number of PS||PS/km2||Area (km2)||Number of PS||PS/km2|
Figure 10 shows the geomechanical classes and the deformation rates on the ore bodies corresponding to the N5W and N5E mines. For both mines, the more relevant ground displacements are related to the poor and very poor geomechanical classes. The relationships between PS and geomechanical classes shown in Table 2 for the mine N5W reveal that the highest amount of detected PS was associated with the largest areas, corresponding to the very poor geomechanical class and reaching displacement rate value up to on the cut slope. This displacement was validated with prism measurements, as discussed in the validation section. For the mine N5E, Table 2 reveals that the poor geomechanical class has the highest area and number of PS, with the high density for all classes, probably due to the low activities in this mine.
Relationship between PS and geomechanical classes for N5W and N5E mines.
|Geomechanical class||Area (km2)||Number of PS||PS/km2||Area (km2)||Number of PS||PS/km2|
In an attempt to validate the results of the proposed approach with field information, a set of measurement points using total station/prisms was used for the validation purpose. To make representative comparisons, the vertical prism measurement data were converted to the LoS of the TSX-1 satellite. The conversion was performed by multiplying the prism measurements by the cosine of the incidence angle (). With the purpose of ensuring that the values obtained by both interferometric and prism measurement refer to the same place, only the closest data pairs were selected. The available field data are related to six sites located along cut slopes of the NW-1 waste pile. These slopes are constituted of highly weathered metavolcanic rocks.26 Figure 11 shows the PSI displacement map, the prism positions, and the related displacement measurements of both techniques. It can be noted that for the period of June to September of 2012 there is an agreement among the different measurements techniques for the six prisms, even considering the prism measurements with a high variability (probably due to variations of total station position of human operation during the measurements), showing any distinguishable deformation trend, which indicates that the slope cuts were stable during the time span. This result was also confirmed by the local geotechnical team of Vale S.A. mining company.
During the period of TSX-1 acquisitions, the slopes of the N5W mine were monitored with total station/prisms. The available prism measurements from April 24, 2012, up to September 28, 2012, were compared with PSI results at the acquisition date from May 3, 2012, up to October 4, 2012, every 11 days, as shown in Fig. 12. The choice of the PS was based on the nearest neighbor of the measured point belonging to the same slope; a weighted average instead could use the PS belonging to different slopes, contaminating the slope displacement itself. The distances from the prism position to the nearest neighbor PS were 1.9, 4.1, and 3.4 m for the prism 1, prism 2, and prism3, respectively.
Topographic measurements by total station/prisms have a high variability [Fig. 12(b)], probably due to the positioning of the total station at each measurement. The computed error between the PSI and prism measurements provided the following results: mean difference equal to 0.085, 1.19 and 1.29 cm for the prisms 1, 2, and 3, respectively, and RMSE of 0.25, 2.29, and 1.95 cm for the prisms 1, 2, and 3, respectively. Even considering these errors, there is a good agreement with PS displacements in terms of trend for the three prism locations.
The proposed methodology took advantage of the previous knowledge and removal of the two more markedly phase components (ground displacement and topographic error) provided by DTS processing; thus, the atmospheric residual phase, orbit refinement, phase noise, and the remaining phase related to the ground displacement and topographic error could be determined more efficiently with the PSI analysis, based on the IPTA technique, making the processing of large areas of extent easier.
The combination of DTS and PSI techniques together with high-resolution TerraSAR-X data, acquired at relatively short intervals (11 day) and covering the time span from March 2012 to April 2013, allowed the detection of linear and nonlinear deformations on the overall mining area. The main affected sectors were the waste piles, reaching accumulated values up to −78 cm, as with the case of the NW-1 waste pile (Fig. 1) (letter A, Fig. 7). Another sector affected by ground displacement was the cut slopes of N5W (letter D in Fig. 7), reaching accumulated values up to . Deep excavations in rock masses of low geomechanical quality coupled with blasting practices and heavy precipitation contributed to this ground displacement. This work presented a methodology that can be applied to monitoring linear and nonlinear ground displacements in a large extent mining area, providing spatial coverage and useful information about ground movement for mining planning and risk assessment.
This investigation was carried out under the scope of an FAPESP-Vale-INPE project (Process FAPESP # 2010/51267-9). The Coordination of Improvement of Higher Education Personnel (CAPES) and the National Council for Scientific and Technological Development (CNPq) are also recognized for a donation received by the first, third, and fifth authors, respectively, during the investigation. The authors would like to thank Vale S.A. for providing access to geological and geomechanical information. Finally, the authors are particularly grateful to the Geotechnical Vale team in Carajás for the support during the field work campaign.
Guilherme Gregório Silva is an environmental engineer. He graduated in environmental engineering from the University of Vale do Paraíba, Brazil, in 2010, and his MS degree in remote sensing from INPE, in 2017. His research interests include topography, cartography, geodesy, environment, remote sensing, GIS, SAR classification, and interferometry.
José Claudio Mura received his BS degree in electric engineering from the University of São Paulo (USP) at São Carlos School of Engineering, Brazil, in 1978, his MS degree in electronics and telecommunications from the Aeronautics Institute of Technology, Brazil, in 1985, and his PhD in applied computer science from INPE, in 2000. He has been with INPE since 1979 and his current research interests include SAR image processing, SAR polarimetry, and interferometry.
Waldir Renato Paradella is a geologist with a BS degree in geology from the USP, Brazil, in 1973, MS degree in remote sensing from INPE, in 1976, and PhD in geology from the USP, in 1983. Since 1974, he has been with INPE and currently is a senior researcher. In 1988 and 1995, he was with CCRS, Canada, as a postdoctoral fellow and visiting scientist. His main research interests include geoscience applications in tropical environments with SAR.
Fabio Furlan Gama graduated in electronic engineering from Vale do Paraiba University, São José dos Campos, Brazil, in 1986, and received his MS and PhD degrees in remote sensing from INPE, in 1996 and 2007, respectively. He currently works at INPE as a senior researcher in the Image Processing Department in development and applications of SAR polarimetry and interferometry. Currently, he is working in estimation of surface deformation using advanced differential interferometric techniques.