Land Subsidence Detection in Tan My-Thuong Tan Open Pit Mine and Surrounding Areas by Time Series of Sentinel-1 Images

Open-pit or underground mining both causes environmental impacts such as air, soil, water pollution, etc., especially causing land subsidence of mines and surrounding areas. Research on mining subsidence is often carried out by field survey, the advantage of this method is high accuracy, but it is usually applied in a small scale. Recently, with the development of radar technology, there have been many studies applying this Radar Interferometry technique to determine surface subsidence over a wide range with a few millimeters accuracy. In this paper, 24 Sentinel-1 images were used as input materials, using the Permanent Scatter Interferometry (PSInSAR) method to determine the land subsidence of the Tan My-Thuong Tan quarries and surrounding areas. The results were compared with the average annual subsidence of 20 surveying points using GNSS technology from 1/2018 to 3/2020. The correlation coefficient of annual average land subsidence of the two methods is 0.83, showing the feasibility of applying InSAR Sentinel-1 data processed by PSInSAR method for determining mine surface deformation and surrounding area.


Introduction
Land subsidence due to resource exploitation such as groundwater and mines is common in many areas of the world such as in Texas (Bawden et  The study of land subsidence by satellite remote sensing has been studied for over 20 years (Chaussard et al. 2014;Chen et al. 2000), in which active radar technology has become a useful tool for monitoring land subsidence.
Differential synthetic aperture radar interferometry (DIn-SAR) was first applied to Seasat satellite imagery to study small elevation changes over an area of 50 km2 in the Imperial Valley, California, USA (Graham 1974). The DInSAR method uses at least two images acquired at two different times of the same location before and after a change occurs to detect the displacement by measuring the different phase of the two times of image acquisition. However, this method has many limitations because it does not eliminate some sources of errors and disturbances, such as atmospheric effects, orbital errors, etc. (Anh et al. 2007).To overcome this restriction, Ferretti proposed the Permanent Scatter SAR Interferometry (PSInSAR) method (Ferretti et al. 2001). The basic principle of this method is based on the use of a series of SAR images in the same area to extract a number of permanent scatter points for determining the land deformations. The PSInSAR method is developed and widely used, which has brought quite good results with accuracy for land deformation detections up to a few millimetres. A typical application such as (López-Quiroz et al. 2009) has successfully applied 38 images of Envisat ASAR to determine subsidence for Mexico City. In Asia, Liu et al. used 26 ERS1 / 2 images to calculate subsidence in Shanghai (Liu et al. 2008). In Indonesia, radar technology was applied for the first time in the land subsidence study of the Jakarta city in 2001 by Hirose's research (Hirose et al. 2001). In this study, 17 SAR images of JERS-1 from February 1993 to September 1998 were used to create 41 interferometric pairs with baselines smaller than 1000 m. Research has shown that during the period of 1993-1995, Jakarta was subsided 10 cm and from 1995-1998 subsidence was 6 cm. In Vietnam, researches to determine urban subsidence is mainly concentrated in Hanoi and Ho Chi Minh City. Anh  with the data to be processed. With Xinpeng Diao research, the authors used advantage DInSAR method with probability integration to detect the settlement on a large scale (Diao et al. 2016), while (Ma et al. 2016) used series of images with SBAS method to determine land subsidence of Bu' ertai Mine, Shendong Coalfield, China. Through these analyzes, it is easy to see the effectiveness of using radar images for monitoring land subsidence in general and mines in particular. Our study focuses on ground subsidence in open-pit quarry mines and surrounding areas by PSInSAR method with Sentinel-1B series from 2018 to 2020. The reason for choosing this type of image is because it is a free, constantly updated image that makes it a great source of data for monitoring mine surfaces. Besides, external measurement benchmarks using GNSS (Global Navigation Satellite System) around the mine were also conducted at the same time of image acquisitions to validate the results for land subsidence of mines and surrounding areas by Sentinel-1 images.

Study area
Binh Duong is a province with many different types of land, suitable for industrial tree development, civil construction and industrial development. In addition, Binh Duong has many special mineral resources such as magmatic rocks, sediments and some weathered rocks. These materials provide the province's traditional industries such as ceramics, construction materials and mining.
The study area is Thuong Tan -Tan My construction quarries located in Thuong Tan commune, Bac Tan Uyen District, Binh Duong province with coordinates ranging from 11o01' to 11° 04' north latitude, and 106° 51' to 106° 54' longitude ( Fig. 2). In the area, there are 17 mines being exploited at different levels, and the deepest is the Chang Tan III and Thuong Tan IV mines currently exploiting at -90 m. According to the provincial mineral exploitation plan, the mines will be assessed the mining capacity down to -150 m. The topography of the research area has an average height of 7 to 40 meters, lowering from north to south, the central area is Lo O Mountain, with an absolute height of 54 meters.

Research methodology
PSInSAR method using a series of multi-temporal images are based on the principle of the DInSAR method. Assuming there is 1 point P on the ground and there are two images that acquired at different times, the phase difference between the two acquisition images can be extracted to determine the land deformations ( Fig. 1). (1) where: M is the first image position (Master), S is the second image position (Slave) B is the distance of the baseline MP is the distance from the satellite at the time of the first image acquisition to the monitoring position P; SP is the distance from the satellite at the second image acquisition time to the monitoring position P; ΔφInt is the interferometric phase determined by the phase difference between two images φM is the measured phase of the first image; φS is the measured phase of the second image; φscatt-M is the phase change generated during the interaction between radar waves of the first acquisition M and target P; φscatt-S is the phase change generated during the interaction between radar waves of the second acquisition S and target P; The above interferometric phase includes the phases related to the terrain element and the deformations together with the noise of the atmosphere, thus to determine the land deformations, it is necessary to remove the terrain phase factor and phase noise. The formula for detecting land deformations is expressed in (2): where: φTopo_res is the residual topographic error (RTE) component; φAtm is the atmospheric phase component at the time of acquisition of each image; φOrb is the phase component due to the orbital errors of each image (errors that affect the position of M and S in (Fig. 1). φnoise is the phase noise; k is an integer value called phase ambiguity, is a result of the wrapped phase of ΔφD-Int, that means the DInSAR phases are bounded in the range (-π, π]; The goal of DInSAR technique is to derive φDispl from ΔφD-Int. This means separating φDispl from the other phase components of Equation (2). An essential condition for doing this separation is to analyze pixels with small φnoise, which often involve to the objects that have strong and constant scattering over time ( Persistent scatteres -PS) and objects also have constant scattering over time, but they come from small scattering objects (Distributed scatteres-DS). The biggest limitation of the DInSAR method is the decrease in correlation as the temporal baseline increases and the φnoise due to atmospheric influences.
PSInSAR method represents an improved method from DInSAR, which uses multiple SAR images acquired in the same area. This method provides suitable data analysis and processing procedures to separate φDispl from different phase components expressed in Equation (2) as well as eliminate the influence of the atmospheric phase screen (APS).

Data and image processing 4.1 Data
The data used is Sentinel-1 images, C band (5.6 cm wavelength). Sentinel-1 satellite operates in four modes (Stripmap (SM), Interferometric Wide swath (IW), Extra-Wide swath (EW), Wave mode (WV)) with different resolutions and has two generations: Sentinel-1A (launched in April 2014) and Sentinel 1B (launched in April 2016). For land surface subsidence studies using the PSInSAR method, the image must be a single look complex (SLC) image, therefore, the selected image mode is IW. Images of the study area are directly downloaded at the Alaska Satellite Facilities (ASF) website of NASA (https://search.asf.alaska.edu/).
In the study area, 24 Sentinel-1B dual-polarization (VV + VH) images with the descending orbit, at the Path of 18 and the Row of 554 have been downloaded. The data is processed for separating the VV polarization. Basic information of images is shown in Table 1. Research area and Sentinel-1B image frame (master image) are shown in Fig. 2. Where the red frame is Sentinel-1B dual-polarized (VV + VH), the small yellow box over the IW1 scanning track is the study area.

Image processing
To process Sentinel-1 images, 2 software: SNAP and StaMPS were used. The process of image processing by PSInSAR method with the two software mentioned above consists of two independent tasks: (i) processing DInSAR for master image and preparing slave images by ESA SNAP; and (ii) PSInSAR processing with StaMPS. Fig.3 depicts a PSInSAR image processing flow chart using SNAP and STaMPS.

Prepare the master image
Firstly, the master image is selected from the download data set, then it is imported into SNAP and selected the sub-swath containing the study area and accurate for the orbit of the Sentinel-1 image by Graph Builder function in SNAP which ability to run the command automatically. This step is important as it will help optimize the time and resources for the rest of the processing (Delgado Blasco et al. 2019). For Tan My -Thuong Tan study areas, IW1 sub-swath was selected (Fig.3).

Slaves preparation
In this step, Sentinel-1 Single Look Complex (SLC) data is sorted by image acquisition dates, and also checks and reduces the name of the original image file. Starting from this processing step to enable processing in batch mode, a Graph Processing Tool (GPT) is used, code for running the processing in xml format is available on Github and run on Python 2.7 (Foumelis et al. 2018).

Image splitting and orbit update
The slave images are also separated based on the part and the polarization, which is determined in the master image preparation step. During further processing, bursts between master and slaves will automatically be extracted according to the defined AOI region. Then the orbital information of slave images are updated by automatically downloading the correct orbital vectors at the ESA Web site (https://qc.senti-nel1.eo.esa.int).

Co-registration and interferogram formation
This step requires a lot of calculations, and it conducts a co-registration of the master image and each slave image that has been prepared in the previous step together. Next, the interferogram for each pair of images was generated and removed the flat earth phase (the phase related to ellipsoid). Debursting for both the SLCs and the differential interferogram is necessary for removing the horizontal stripes. The terrain phase is simulated and removed using the 3-second SRTM Digital Terrain Model (DTM), which is automatically downloaded by SNAP. In this step, the research area that has reserved before is subsetted.

Export to StaMPS format
This is the last step of the DInSAR processing chain on SNAP. Exporting this data will create folders containing single images of SLC views of all image files, folders containing the interferometric pairs of the master image and the slave images as well as the folder of the coordinates of the master image and the study area digital elevation model (DEM).

Ingestion into StaMPS
The following step is the ingestion of SNAP exports into StaMPS using a specific script, called mt_prep_snap (Foumelis et al. 2018). The StaMPS PSInSAR processing sequence is then run from steps 1 to 7 (Hooper et al. 2010). Below, we explain a little more about some key steps in the STaMPS image processing chain.

PS Selection
This step makes a selection of PS based on probability, by comparison to results for data with random phase. The main part of the time series calculations are done in Matlab (Hooper et al. 2010).

Phase unwrapping
The interferometric phase when receiving is wrap phase, and it varies from -π to π, so in order to get the real phase of the terrain, we have to do phase unwrapping. The phase unwrapping technique is the most difficult step and also the decisive step to the accuracy of the results of determining the terrain deformations made by InSAR method. StaMPS has the phase unwrapping methods as the "Minimum Cost flow 2D (MCF)" algorithm, or the MCF 3D method (Hooper et al. 2010). The MCF 3D method was chosen to the phase unwrapping for our data set because it proved highly accurate (Hooper et al. 2010).

Subsidence calculation over time and eliminate atmospheric influences
The heterogeneity of the atmosphere (ionosphere and troposphere) and its variation over time and space, cause changes in the signal speed on the line between the antenna and the terrain surface, and that is the cause that directly affects the interferometric phase value and is called the Atmospheric Phase Screen (APS). In this step, we will have to remove the APS first and calculate the speed of land subsidence.
The speed of subsidence is calculated from a series of subsidences over time of different times. Suppose subsidence of each target is d = [d1, d2,. . ., dn] (n is the number of the image acquisitions) and the baseline of the images in time is T = [T1, T2,. . . Tn], then the weight is used to calculate the rate of subsidence, provided that the mean square error of the interferometric phase will be taken as the weight of the settlement. The following formula is used to determine the rate of subsidence over time (Jiang et al. 2016). where: υ is the subsidence velocity, T represents the corresponding time interval between image acquisitions d is a time series deformation, P is a weight matrix and is defined as:

Prepare the master image
Firstly, the master image is selected from the download data set, then it is imported into SNAP and selected the subswath containing the study area and accurate for the orbit of the Sentinel-1 image by Graph Builder function in SNAP which ability to run the command automatically. This step is important as it will help optimize the time and resources for the rest of the processing (Delgado Blasco et al. 2019). For Tan My -Thuong Tan study areas, IW1 sub-swath was selected (Fig.3).

Slaves preparation
In this step, Sentinel-1 Single Look Complex (SLC) data is sorted by image acquisition dates, and also checks and reduces the name of the original image file.
Starting from this processing step to enable processing in batch mode, a Graph Processing Tool (GPT) is used, code for running the processing in xml format is available on Github and run on Python 2.7 (Foumelis et al. 2018).

Image splitting and orbit update
The slave images are also separated based on the part and the polarization, which is determined in the master image preparation step. During further processing, bursts between master and slaves will automatically be extracted according to the defined AOI region. Then the orbital information of slave images are updated by automatically downloading the correct orbital vectors at the ESA Web site (https://qc.senti-nel1.eo.esa.int).

Co-registration and interferogram formation
This step requires a lot of calculations, and it conducts a co-registration of the master image and each slave image that has been prepared in the previous step together. Next, the interferogram for each pair of images was generated and removed the flat earth phase (the phase related to ellipsoid). Debursting for both the SLCs and the differential interferogram is necessary for removing the horizontal stripes. The terrain phase is simulated and removed using the 3-second SRTM Digital Terrain Model (DTM), which is automatically downloaded by SNAP. In this step, the research area that has reserved before is subsetted.

Export to StaMPS format
This is the last step of the DInSAR processing chain on SNAP. Exporting this data will create folders containing single images of SLC views of all image files, folders containing the interferometric pairs of the master image and the slave images as well as the folder of the coordinates of the master image and the study area digital elevation model (DEM).

Ingestion into StaMPS
The following step is the ingestion of SNAP exports into StaMPS using a specific script, called mt_prep_snap (Foumelis et al. 2018). The StaMPS PSInSAR processing sequence is then run from steps 1 to 7 (Hooper et al. 2010). Below, we explain a little more about some key steps in the STaMPS image processing chain.

PS Selection
This step makes a selection of PS based on probability, by comparison to results for data with random phase.

Phase unwrapping
The interferometric phase when receiving is wrap phase, and it varies from -π to π, so in order to get the real phase of the terrain, we have to do phase unwrapping. The phase unwrapping technique is the most difficult step and also the decisive step to the accuracy of the results of determining the terrain deformations made by InSAR method. StaMPS has the phase unwrapping methods as the "Minimum Cost flow 2D (MCF)" algorithm, or the MCF 3D method (Hooper et al. 2010). The MCF 3D method was chosen to the phase unwrapping for our data set because it proved highly accurate (Hooper et al. 2010).

Subsidence calculation over time and eliminate atmospheric influences
The heterogeneity of the atmosphere (ionosphere and troposphere) and its variation over time and space, cause changes in the signal speed on the line between the antenna and the terrain surface, and that is the cause that directly affects the interferometric phase value and is called the Atmospheric Phase Screen (APS). In this step, we will have to remove the APS first and calculate the speed of land subsidence.
The speed of subsidence is calculated from a series of subsidences over time of different times. Suppose subsidence of each target is d = [d1, d2,. . ., dn] (n is the number of the image acquisitions) and the baseline of the images in time is T = [T1, T2,. . . Tn], then the weight is used to calculate the rate of subsidence, provided that the mean square error of the interferometric phase will be taken as the weight of the settlement. The following formula is used to determine the rate of subsidence over time (Jiang et al. 2016). (3) where: υ is the subsidence velocity, T represents the corresponding time interval between image acquisitions d is a time series deformation, P is a weight matrix and is defined as: (4)

Results of determining the land subsidence from the Sentinel-1 time-series images
Results of Sentinel-1B images processing for Tan My-Thuong Tan area after exporting from STaMPS is imported to ArrcGIS10.2 (Fig.4). In the research area, there are many subsidence points, and the average deformation value is less than 15 mm/1 year. The number of points with the high rate of subsidence is about -15 mm/year while the majority of subsidence is smaller than -10 mm/year. However, the sites around the mine are located on roads or houses that are our interest because these sites will directly affect people's lives.

Validation of subsidence results
In order to validate the results of determining land subsidence from Sentinel-1 satellite images, 20 monitoring points were designed and measured over the entire area at locations that can reflect the deformations of the Tan My -Thuong Tan area. Because of its large area, the measurement method chosen was the GNSS high accuracy. The monitoring benchmarks were connected with one benchmark of 1st state grade, and this benchmark was measured and compared to the other two elevation benchmarks to assess the stability before each monitoring cycle, which was not displacement in both measurement cycles in 2018 and 2020. The static measurement method was chosen with the measurement devices were five receivers of GNSS CHC X91B. The number of measurements  Table 2.

Discussions a. At mining areas
The subsidence points from Sentinel-1 images are quite concentrated in the mine area. However, because it is an open-pit mine, the surface of the mine topography is constantly changing with the exploitation process, so that the mining pits have not been surveyed in the field by GNSS.
Although there is not any GNSS measurement data to validate subsidence at the mine pits, from the comparison results in Table 2 which is the subsidence survey sites by GNSS and the PS points located near the quarry or on roads close to the open pits, it can be seen that the subsidence values from the images are quite similar to the result from the GNSS measurement, with accuracy ± 1.5 mm. Therefore, the results of subsidence at the open pit sites are considered accurate. With the determined values of subsidence in the mine area, we can note some locations with large values of subsidence, such as at locations A, B and C (Fig.6) which are concentrated on the pits and the transport routes.
For A position, the maximum value of subsidence is about -13 mm/year with red PS points, orange PS points with subsidence values less than -10 mm/year. At B position, the number of distribution subsidence points are more than the A position, but the subsidence values are not as high as those in A position (see the chart in Fig.6). Subsidence values are mostly less than -5 mm/year, some are less than -10 mm/year, and only 4 PS points are less than -15 mm/year, and some are positive. At C position, the average of subsidence is about -10mm/year, there are some PS points with subsidence close to -20mm/year. The majority of subsidence ranges from -5 mm to -10 mm/year.
The reason for the occurrence of PS points with both positive and negative values in the mining pits is partly due to the drilling and blasting operations, and the work of heavy trucks are also a cause of the positive deformation.  Table 2 above. The reason for the selection of points within this 20 m radius is due to the scattering properties of the radar image when the Radar wave reaches the surface of the object, the reflected scattering rays may not return immediately, but it may be corner backscattering or volume scattering before returning, thus the checkpoints shown on the image may be slightly shifted from its actual position.
Based on Table 2, the GNSS measurements have 2 years interval while the measurements from 24 Sentinel-1 image processing are averaged over the period from January 2018 to March 2020. So the values of subsidence measurements using GNSS will have to be divided by an average of 2 years. The correlation of these 20 points was evaluated for the purpose of a preliminary assessment of how the subsidence points measured by the PSInSAR method achieved compared to the high precision GNSS measurement method. Below is the distribution of ground settlement values of the 20 locations shown in Table 2.
Based on the chart in Fig. 7, it is to recognize that the settlement values made from PSInSAR method tend to be higher than the settlement values measured from GNSS high accuracy method; however, the correlation is quite good with R2 reached 0.83. This can also be explained easily because the number of survey points is not much, and in addition, the GNSS and PS points from the image are not completely overlapping. Although the number of checkpoints is not much (20 points), the measured values also reflect the subsidence situation around the Tan My and Thuong Tan open-pit mines in the period from January 2018 to March 2020

Conclusions
With a data set of 24 Sentinel-1B images collected between January 2018 and March 2020, the PSInSAR method was applied with the combination of ESA SNAP and StaMPS software. The results of determining the land subsidence of Tan My -Thuong Tan open pit and surround area indicate the following: The PSInSAR method is suitable for mining areas with few trees and as many constructions as this study area. Sentinel-1 image with a resolution of 3.5 m in range and 22 m in azimuth, large cover and high frequency of image acquisitions (12 days), is very appropriate for studies related to land subsidence.
This is the first time, the issue of land subsidence in mining areas in Vietnam has been investigated by using a combination of radar interferometric method and GNSS method. With a large number of images, the rate of land subsidence was evaluated within a period of approximately 2 years and 3 months. In the quarry area, there were many subsidence PS points, but these subsidence points were mostly located on the road of pits or the slope of the mining pits. The biggest land subsidence rate at mines was smaller than 15mm / year.
In the area around the quarry such as roads, the edge of pits or houses near the mines have been surveyed by GNSS. PSInSAR subsidence points were compared to GNSS measurement points from January 2018 to March 2020. However, the number of GNSS measurement points and PS points from radar images did not completely coincide, and we used 20 locations around the mine area where PS subsidence points coincided with those measured by GNSS in a radius of 20 m. Correlation between these two data types reached 0.83, proving the ability to determine land subsidence by the Sentinel-1 image sequence.