Pyrenees deformation monitoring using Sentinel-1 data and the Persistent Scatterer Interferometry technique

This study focuses on deformation mapping and monitoring using Sentinel-1 radar data and the DInSAR (Differential Interferometric Synthetic Aperture Radar) technique. In particular, a Persistent Scatterer Interferometry technique was used over 15000 𝑘𝑚 2 of the Pyrenees, located in Spain, Andorra and France. The main goal is to monitor deformations over a vegetated and mountainous region by exploiting the wide-area coverage, the high coherence and temporal sampling provided by the Sentinel-1 satellites. All possible interferograms were generated using 150 images covering five years. The velocity map of the entire region is presented considering the characteristics of the study area, including vegetation and severe steep mountains. Two areas of deformation are shown, which are characterized by velocity values ranging between -20 to -40 mm/year.


Introduction
Remote sensing has an extensive capability to monitor the nature of Earth and the human-made activities. Spatial and temporal data of different satellites enable researchers to investigate changes such as land deformation [1]. Land deformation may cause ground surface ruptures, increasing flood risk, infrastructure damage, and adverse socioeconomic impact in the community [2]. Interferometric Synthetic Aperture Radar (InSAR) is a powerful technique that precisely measures ground surface deformation over wide areas. Spaceborne differential interferometric SAR (DInSAR) is one of the most applicable methods to monitor land deformation. Satellite-based SAR systems provide very high spatial and high temporal resolution (4-6 days nowadays) over a large area in all weather conditions, day and night [3]. Persistent Scatterer Interferometry (PSI) is a powerful DInSAR technique able to measure sub-metric topography and millimetric deformation [1,4]. PSI represents a specific class of DInSAR techniques exploiting the information contained in the radar phase of at least two complexes SAR images acquired in different times over the same area, which are used to form an interferometric pair [1]. First presentation of DInSAR was provided at 1989 based on L-band Seasat data [5]. A general review of SAR interferometry and probable applications have been discussed in Bamler and Hartl, 1998 [6] and Rosen et al., 2000 [7]. The first PSI technique is the Permanent Scatterers approach proposed by Ferretti et al., 2000 [8] and Ferretti et al., 2001 [9]. Several other contributions then followed these pioneering works. The Small Baseline Subset (SBAS) technique [10] is one of the most extensively used. Based on multi look imagery, Mora et al., 2003 [11] and Schmidt and Bürgmann, 2003 [12] investigated similar approaches to SBAS. Crosetto et al. (2005) described a simplified PSI approach based on stepwise linear deformation functions and least-squares adjustment [13]. To exploit the relative properties of neighbouring pairs of PSs, Costantini et al. (2008Costantini et al. ( , 2014 [14] and [15], proposed another PSI method. Crosetto et al. (2016) provided an extensive and complete review of persistent scatterer interferometry [1], including the main PSI algorithms, several technical aspects of the PSI techniques, the validation of the PSI results, the main PSI applications, and main open PSI problems. Finally, [16] and [17] investigate landslides in reservoirs of southern Spain and over wide areas with geological threats of Italy, based on the developed technique of [4].
In this study, an approach of Earth's surface deformation using SAR images of a vast vegetated and mountainous region located over the Pyrenees is presented. The organization of the paper is as follows. Section 2 first provides information on the Sentinel-1 data; then, it describes the situation of the study area. An overview of the used PSI method is explained in Section 3. The velocity maps are shown in Section 4. Finally, Section 5 presents a conclusion of the study.

Study area
The study area includes 15000 2 of the Pyrenees, located in Spain (Catalonia), Andorra, and France ( Fig. 1). This area contains high mountains with steep terrains, rising higher than 3,400m in elevation. The area is mostly snowy during winter. Moreover, the land cover of the region mainly includes vegetation, forest, small urban areas, and high mountains. Landslides and avalanches are the main natural hazards. The climate varies significantly with elevation. The valleys have a climate that is similar to the temperate climate of Andorra's neighbours, but because of the higher elevation, in the mountains winters tend to be more severe, the humidity lower, and summers slightly cooler. The average annual precipitation varies across the region, increasing with elevation and from south to north. The driest months tend to be January and February, and the wettest, May, June and November. During the summer months, there are exceptionally few rainy days, but the rainfall can be very weighty because it is associated with thunderstorms.

Dataset
In Table 1, the main characteristics of the used dataset are described. The dataset consists of 150 Sentinel-1 Interferometric Wide Swath images covering around five years, with the first acquisition time in March 2015 and the last acquisition time in November 2019. In this study, images from Sentinel-1A and B satellite have been used; thus, the minimum temporal sampling is 6 days. 1500 interferograms have been generated to derive the deformation maps.  A single IWS frame of Sentinel-1 covers entirely the study area. From this frame, which covers approximately 240 by 180 km, twelve bursts coming from two swaths were selected and extracted. The interferometric processing, which is described in the next section, was performed on the twelve bursts. In order to compute the differential interferometric phase, a SRTM Digital Elevation Model with 30 m cell resolution was used to simulate and to remove the topographic phase component.

Method
PSI [1] is an opportunistic deformation measurement method i.e., it can measure deformation only over the available Persistent Scatterers (PS). The PS density is usually low in vegetated, forested, and low-reflectivity areas (e.g., very smooth surfaces), and in steep terrains facing the radar sensors. Snow coverage, construction works, street repavement, etc. can cause the complete or partial loss of PSs. By contrast, PSs are usually abundant in buildings, monuments, antennas, poles, conducts, exposed rocks, or outcrops, among others. In this research, the deformation maps have been generated using the Persistent Scatterer Interferometry chain of the Geomatics (PSIG) Division of CTTC, which is described in [18].
The main steps of the processing used in this work are briefly described in the following. It is worth noting that this represents a simplified version of the PSIG chain.
1) The interferogram network is generated using the Sentinel-1 images. Since a fundamental aspect of the PSIG chain is the redundancy of the network of interferograms and images, all the possible interferogram pairs are usually generated.
2) For those areas where the winter snow coverage is important, it is necessary to remove from the stack of interferograms the affected interferograms. The selection of the interferogram network is usually made by evaluating the coherence images.
3) Only a portion of the millions of pixels available in a single frame is exploited for deformation measurement. The main goal of point selection is to find a reasonable compromise between the quality of the selected points (little affected by noise) and a proper spatial coverage. The point selection can be based on the amplitude dispersion [8] or the coherence. 4) Over the selected pixels and using the redundant network of wrapped interferograms, the linear annual deformation velocity and the residual topographic error are estimated. This operation, which is based on the periodogram, is first performed on a set of first order points. The main goal of this operation is to cover homogeneously most of the area of interest.
5) The first order processing is followed by a second order analysis, whose main objective is to densify the first order network of points.
6) The final output is a velocity map including, for each point, the estimated deformation velocity in the satellite Line of Sight direction. Average of accumulated values has been used to mitigate the atmospheric phase screen effects.

Results and Discussion
Starting from the generated interferograms and using the analysis of the coherence, which was mainly focused on decorrelation due to snow cover, 8452 interferograms were discarded. Using the remaining 814 interferograms, the linear deformation velocity was calculated over the selected points. Note that the linear deformation velocity is estimated along the satellite Line of Sight (LoS). A general overview of the velocity map, covering an area of around 15000 2 , is shown in Figure 2. The velocity map includes 2904322 measured points. The green colour, with velocity values ranging from -3.5 to +3.5 mm/year, represent areas were no movement has been detected.
In the area shown in Figure 2 there are hundreds of localized deformation phenomena. This shows the great potential of this technique over wide areas. In this work, we only mention two detected slope movements, whose location is indicated by a red box in Figure 2. These two areas are affected by deformation. Figure 3 shows a zoom of the two areas, both are located in the Ariège department (Occitanie region, France): one in thevalley of Artiès, and the other one close to the lake of Soulcem. To assess the reliability of the results, specialists of ground-based land deformation and geology have been analyzing the obtained results.
The area shown in Figure 3B (Arties) is globally stable. An exception is given by a group of yellow and red pixels. The interpretation of this group is that there is clearly a landslide phenomenon, which affects one road. The most active part of the landslide has displacements above the 15 mm/year, reaching in some areas 40 mm/year. This is a key information for those that manage natural hazard risk.
Another interesting example is given by the Lake Soulcem, where a set of deformation phenomena occur in a flank that potentially affects the artificial lake. The main deformation phenomenon interests an important area with the most active zone reaching a velocitiy of 20 mm/year.

Conclusion
PSI is a remote sensing technique able to measure millimetric displacements of the Earth's surface from space. This paper aims at explaining the implemented methodology to generate deformation maps over broad areas, using the DInSAR technique and Sentinel-1 SAR data. The outputs of the methodology are deformation maps in terms of velocity maps. The application and the results of the methodology over the Pyrenees were presented and discussed. The derived deformation map covers an area of approximately 15000 km 2 . It includes hundreds of deformation phenomena: this shows the great potential of the technique. This potential can be fully exploited considering the key characteristics of the PSI data, and integrating the PSI information with other types of information, e.g. geology maps, digital terrain models, etc.