LANDSLIDES MONITORING WITH TIME SERIES OF SENTINEL-1 IMAGERY IN YEN

1Dept. of Geomatics and Land Adminnistration, HUMG, HaNoi University of Mining and Geology, No 18 Vien Street, Bac Tu Liem, Hanoi, Vietnam tranvananh@humg.edu.vn 2Dept. of Information Technology, HUNRE, HaNoi University of Natural Resources and Environment, No 41 A Phu Dien Street, Bac Tu Liem district, Hanoi, Vietnam Txquang@hunre.edu.vn 3Institute of Geological Sciences, Vietnam Academy of Science and Technology, No. 84 Chua Lang Street, Dong Da, Hanoi, Vietnam nguyenducanh237@gmail.com 4Department of Civil and Environmental Engineering (DICA) Politecnico di Milano, Piazza Leonardo da Vinci 32, Milan, Italy(laura.longoni,vasil.yordanov)@polimi.it 5 Vasil Levski National Military University, Veliko Tarnovo, Bulgaria


INTRODUCTION
Vietnam is one of the countries in the region that is frequently affected by landslides due to tropical monsoon climate and three-fourths of Vietnam land area is hilly. In the context of global climate change happening quite acute, landslides are becoming more dangerous and more severe. According to research by (Quoc Hung et al., 2017) almost every year in Vietnam during the rainy season landslides occur, causing great damage to people and properties.
From the early years of the twenty-first century, scientists around the world have studied the problem of landslides and published many valuable papers on this field. Many works have used remote sensing data to detect landslide regions, tectonic destruction zones, etc., In the remote sensing field, the SAR satellite interferometric measurement is a method applied for evaluating changes on the Earth's surface proved very useful. Differential SAR Interferometry (DInSAR) is a method used to detect land deformations by using two or more images of the same location acquired at different times, usually before and after a topographic change occurs (Anh et al., 2007). However, this method has many limitations such as the influence of the atmosphere and some scattering characteristics of objects on the surface. The PSI method proposed by (Ferretti et al., 2001), is based on using a series of multi-temporal SAR images of the same location to extract the permanent scattering points which are used for detecting terrain deformation.
PSInSAR method has been applied to monitor landslides by (Colesanti et al., 2003), (Colesanti and Wasowski, 2006), (Bozzano et al., 2017). In particular, the existence of free SAR images in the archives as well as new acquisitions gives PSInSAR the ability to measure and monitor the terrain changes in the past and present.
There have been many landslide studies in Vietnam, including a study on landslides in Hoa Binh Hydropower (Hung and Tung, 2015). Another research about predicting the risk of landslide in mountainous areas, North Vietnam and solutions for prevention by Dieu Tien Bui (Bui et al., 2017). Recently, the research of (Van Tran et al., 2020) used ALOS PalSAR images with Small Baseline (SBAS) method for detecting landslides in Laocai province. In our research, the PSInSAR method was used for landslide determination in the Van Yen district, Yenbai province of Vietnam using the Sentinel-1 time series images.

METHODOLOGY
The use of radar interferometry (InSAR) technique to monitor the surface of the Earth, including the topographic surfaces and the terrain deformation, has been successfully demonstrated over the past two decades. The traditional interferometric method allows the creation of interferogram images of the phase shifts between two or three images acquired at different times over the same area on the surface. Equation (1) illustrates the phase of land deformation (Zebker et al, 1986).
Assuming that if a DEM of the imaged scene is available, φTopo can be simulated and subtracted from ΔφInt (this is the inverse operation performed in InSAR DEM generation), obtaining the so-called DInSAR phase ΔφD-Int: The goal of DInSAR technique is to derive φDispl from ΔφD_Int. This implies separating φDispl from the other phase components of Eq. (1). An essential condition to accomplish this separation is to analyse pixels characterized by small φNoise, which are typically related to two types of reflectors: those where the response to the radar is dominated by a strong reflecting object and is constant over time (Permanent Scatterer, PS) and those where the response is constant over time but is due to different small scattering objects (Distributed Scatterers, DS). PSI represents a specific class of DInSAR techniques, which exploits multiple SAR images acquired over the same area, and appropriate data processing and analysis procedures to separate φDispl from the other phase components depicted in Eq. (1). For more information on PSInSAR refer to Ferreti's documentation (Ferretti et al., 2001).

Study area
The study area is Van Yen district which is the largest district in Yen Bai province, Vietnam. There is provincial road 151 and the Hanoi -Lao Cai railway line running along the area connecting Tran Yen district with Lao Cai province.
Van Yen has a very rich system of rivers, streams, ponds and lakes. The Red River originates from Yunnan (China), the length flowing through Van Yen is 70 km long. The tributaries of the Red River in the district have up to 40 streams flowing into the Red River. The largest of which are Thia and Hut streams flowing from Van Chan district through the district with a total length of more than 100 km.
The topography of Van Yen is mainly in the form of medium to high mountains, originating from corrasion and denudation, the elevation difference is from 200 m to 1700 m, gradually lower to the southeast. The whole Van Yen district is like "a mountain trough" that extends from the northwest to the southeast, in the middle is the Red River valley. This area has a natural slope from 20 o to 40 o with a network of many streams, narrow steep valleys and is strongly dissected. This is one of the very important factors causing landslides and other disasters in this district ( Figure 1).
The study area has a very complex geological structure, neotectonic activities develop strongly, with the major fault being the deep fault of the Red River which forms a fault system that develops in the northwest -southeast. These fault systems have strongly changed the structure of the bedrock, creating many tectonic rupture zones, which are very favourable conditions for occurring of geological disasters. Through field survey, it was found that Provincial Road 151 running through almost coincides with the faults of the Red River. Combined with the process of sawing the slopes to make roads, many landslides occurred on Provincial Road 151 and that is the cause for greatly affecting the life, economy -society of people in Van Yen district.
As reported by Le Quoc Hung (Quoc Hung et al., 2017), the phenomenon of landslides and geological disasters in this area often occurs in the rainy season from June to October every year. Landslides often occur on talus cliffs or mountain slopes with human impact such as cutting slopes to make roads. At present, landslides happen in many and densely along the road section 151 through Dong An, Chau Que Ha, Dong Cuong, Lam Gang communes. (Figure 2). In addition, there are large-scale landslides on the mountains due to tectonic activities. Figure 2 shows the survey location of landslide sites, last updated in 2013.

Data
In this study, 27 descending Sentinel-1A (Terrain Observation with Progressive Scans SAR)-TOPSAR images, C band and (VV+VH) polarization acquired from 11 th January 2019 to 1st March 2021 were used. Sentinel-1A are free accessibility at Alaska Satellite Facility archives (ASF) with regularly repeated acquisition at a 12-day interval (https://asf.alaska.edu/ ). Details of the Sentinel-1A data used in the study including the date of acquisitions and the perpendicular baseline are shown in Table 1. All pairs of SAR images were checked for concordance by selecting a perpendicular baseline of less than 200m (Tzouvaras et al., 2020).

IMAGE PROCESSING
With the Sentinel-1 dataset, the data processing for PSInSAR method will be divided into two main parts (Foumelis et al., 2018): (I) DInSAR processing for the master image and preparation of the slave images for interferogram generation using SNAP software (ESA SNAP), and (II) PSInSAR processing using StaMPS (Hooper et al., 2012). Figure 3 is a flow chart of PSInSAR data processing using SNAP and StaMPS. Details of these steps will be described in the following sections. Figure 4 shows the temporal and perpendicular baseline of data set, in which the central position is the master image and the surrounding positions are the slave images.

Select the sub-swath for master image
After selecting the main image from the dataset in SNAP software, this image will be extracted to a sub-region containing the study area because Sentinel-1 image has 3 sub-swaths in the vertical direction and 9 bursts in the horizontal direction. This The International Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences, Volume XLVI-4/W2-2021 FOSS4G 2021 -Academic Track, 27 September-2 October 2021, Buenos Aires, Argentina selection is run automatically using the Graph Builder function in SNAP. This step is quite important because it reduces the time and memory for the rest of the processing (Delgado Blasco et al., 2019). For the Van Yen study area, the IW1 scan was selected.

Slave image files preparation
Except for the image selected as master image, the remaining images in the Sentinel-1 data set will be sorted in the order of the date of image acquisition and included as slave images. From this processing step, batch mode processing is possible, the graph processing tool (GPT) is used, the code to run the processing in XML format is available on Github and runs on Python 2.7 (Foumelis et al., 2018).

Spliting the slave images
The slave images are also depolarized as VV is similar to master image. During further processing, master and slave images are automatically subsetted according to the defined the area of interest (AOI). The slave image orbital information are then updated by automatically downloading the correct orbital vectors from the ESA web site (https://qc.sentinel1.eo.esa.int)

Image co-registration and interferometric generation
This step is computationally demanding, since the master image and each slave image, previously prepared, are co-registered. Next, interferogram for each pair of images is generated and the flat earth phase is removed (the phase relative to the ellipsoid). Debursting all the SLC images is necessary for removing the horizontal stripes. The terrain phase is simulated and eliminated using the SRTM Digital Terrain Model (DTM), which is automatically downloaded by SNAP. In this step, the study area will be subset.

Export the interferometric images into STaMPS
Exporting the data will produce three folders: first containing single look complex images (SLC) of all image, second containing the interferometric pairs of the master image and the slave images and the last folder is containing the coordinates of the master image and the digital elevation model (DEM) of the study area.

PS selection
After exporting the data from SNAP, it needs a preparatory step before importing it into STaMPS. This preparation helps for the selection of permanent scattering points (PSs) of the time series images. The threshold for PSs point taken depends on the amplitude dispersion DA that is shown as: ( 2) where ¶A is the standard deviation and mA is the mean of the backscattering intensity. The main part of the PSInSAR calculation in STaMPS is done in Matlab. The selection of PS points of the images are done according to the DA value set from the previous step, then the image pair coherence is calculated and selected the PSs (Hooper et al., 2010).

Phase unwrapping
The interferometric phase of image pair after calculation is the wrap phase which varies from -p to p and is called ambiguity. Phase unwrapping is used to remove ambiguity in order to get real phase. Phase unwrapping is the most challenging step and also it could influence the accuracy of the results. 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 for the phase unwrapping of our dataset because it proved highly accurate (Hooper et al., 2010).

Landslides calculation over time and eliminate atmospheric influence
In this step, spatially-correlated look angle (SCLA) error is computed which is due to almost exclusively to spatiallycorrelated DEM error. The master atmosphere and orbit error (AOE) phase is estimated simultaneously (Hooper et al., 2010).

RESULTS AND DISCUSSION
With the dataset downloaded and using the SNAP and StaMPS software, it was obtained 54666 permanent scattering points (PSs) in the study area. After phase unwrapping and atmospheric influence removed, the displacement value in the light of sight (LOS) was calculated and exported to shape file. This result is superimposed on the topographic map including roads and river system ( Landslides survey points collected in 2013 by Vietnam Institute of Geosciences and Mineral resources (Quoc Hung et al., 2017) were mapped and compared with some positions of landslide points analysed from the radar images (PSs). There were some spots that are quite similar to the surveying results. The areas with dense PS points will be investigated in more detail below are the communes of Lang Thip, Chau Que Ha, Xuan Tam, Phong Du Thuong and the area along provincial road 151.

Lang Thip commune landslide
Lang Thip commune has a rather complicated topography, mountains with an average height between 900 m and 1100 m. The PS points are distributed mainly on steep slopes along inter-commune roads. The time series for Sentinel-1 data from January 11, 2019 to March 1, 2021 ( Figure 6) is not smooth and shows the displacement rate in the LOS direction of PS points in this area of average -15.4 mm/ year. In the region, there are not landslide survey points during this time period except for the survey data from 2013 mentioned above. Based on the field survey points overlaid on the landslide map from SAR analysis, the field survey point with the higher deformation value inside area A (YB120252.L2) was selected and graphed. The graph shows land deformation in the time from January 11, 2019 to March 1, 2021 ( Figure 6).

Chau Que Ha commune landslide
In the Chau Que Ha area, PS points are concentrated in 3 clusters, these points also have medium to small landslide velocty, and the average landslide rate is less than -10mm /y. In this area, the landslide survey points are not very matching, in two clusters of high deformation, there are no survey points. At a survey site inside ellipse B with code YB120583.L2 is matched with PS point. This point has been charted by time series (Figure 7). The variation of the soil deformation values is quite coarse, which ranges from +20mm to -20mm.

Xuan Tam commune landslide
Similar to the Chau Que Ha, clusters of PS points appear quite a lot on slopes with slopes from 15 to 20%, PS points have a landslide rate of less than 10mm/y in the LOS direction ( Figure  8). In this commune, there are four clusters of deformation points clustered in four ellipses. In these four clusters, the field survey points also do not coincide with the PS points but shift with each other. In ellipse C, a point has been taken close to the survey point with code YB120491.L2 The displacement velocity graph is reported in Figure 8c.

Phong Du Thuong commune landslide
Phong Du Thuong is located in the southwest of Van Yen, this area has an average altitude of 200 m to 1500 m. The distribution of PS points here is quite concentrated in an area along the inter-commune road and the streams. The location with many PS points has been highlighted by the ellipse. Unfortunately, there is almost no field survey in this area, except for a single point, but it does not coincide with any PS point.
Because the PS points showing the magnitude of landslides appear to be the highest in Van Yen district, a location within The International Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences, Volume XLVI-4/W2-2021 FOSS4G 2021 -Academic Track, 27 September-2 October 2021, Buenos Aires, Argentina the circle in Figure 9b were selected for charting and the graph is shown in Figure 9c. The graph shows landslide values over time from January 2019 to March 2021 which are +40mm to -40mm.

Na Hau commune landslide
In the Na Hau, deformation values are small below 10mm/year, concentrated in the road connecting Dai Son commune to the centre of Na Hau commune. In the area illustrated by two ellipses A and B, the average landslide rate are less than 10mm/y and 13mm/y at A and B location, respectively. In ellipse B, there is one field survey point with code YB160197.L2 that coincides with some PS points and it is used to report the velocity chart shown in Figure 10c.

Along provincial road 151 landslide
For the area along Provincial Road 151, the terrain is not characterised with high altitudes, with an average elevation of less than 200m. Numerous PS points appear, but all values show small deformation, most of the average land deformation values are less than 5mm/year except for a few higher values which are matched with the field survey points from 2013. They those are located at Chau Que Ha, Lam Giang, Dong An with assigned code numbers YB120336.L2, YB120614.L2, YB120558.L2. From these locations, a point in Lam Giang commune was selected to plot the deformation over time and the results are shown in Figure 11c. When interpreting landslide maps of the study area based on SAR Sentinel-1 data, it is important to note that not all slopes are visible by SAR due to the way of the image acquisition causes shadow and layover. In addition, the exposure of the landslide slopes is also important because the motion of the scattering points in the direction close to the satellite track of the satellite can hardly be corrected by the SAR interferometer (Kiseleva et al., 2014). This explains the fact that there are many locations where the landslide occurred but could not be identified by radar images. In addition, dense vegetation also limits the backscattering of radar waves, so it is difficult for areas with dense forests like Xuan Tam, Chau Que Ha.

CONCLUSION
The use of multi-temporal Radar satellite imagery with PSInSAR method helps to better understand and model the progress of landslides. In addition, based on the time-series of Sentinel-1 data the landslide velocity can be generally calculate at certain times even for very small deformations (mm / year).
However, this method also has some restrictions since they heavily depend on the choice of satellite orbits are suitable to the area. For some locations landslides can be identified in an ascending direction however in a descending direction it is not possible. Therefore, land deformations can be more accurately determined if both ascending and descending images are processed.
For areas covered with dense vegetation there are no permanent scatter points (PS) that can be monitored at various times or another issue could be the size of the PS points not being large enough, making this method hard to apply. To overcome the problem, artificial scattering points should be installed in the study area, which will act as PS points and the reflected signal for the points would be very good.