Open AccessArticle

Automatic Generation of Sentinel-1 Continental Scale DInSAR Deformation Time Series through an Extended P-SBAS Processing Pipeline in a Cloud Computing Environment

Riccardo Lanari

Manuela Bonano

Francesco Casu

Claudio De Luca

Michele Manunta

Mariarosaria Manzo

Giovanni Onorato

and

Ivana Zinno

Institute for Electromagnetic Sensing of the Environment, National Research Council (IREA-CNR), via Diocleziano 328, 80124 Naples, Italy

Author to whom correspondence should be addressed.

Remote Sens. 2020, 12(18), 2961; https://doi.org/10.3390/rs12182961

Submission received: 27 July 2020 / Revised: 7 September 2020 / Accepted: 9 September 2020 / Published: 11 September 2020

(This article belongs to the Special Issue Scaling-Up Deformation Monitoring and Analysis)

Download

Browse Figures

Versions Notes

Abstract

We present in this work an advanced processing pipeline for continental scale differential synthetic aperture radar (DInSAR) deformation time series generation, which is based on the parallel small baseline subset (P-SBAS) approach and on the joint exploitation of Sentinel-1 (S-1) interferometric wide swath (IWS) SAR data, continuous global navigation satellite system (GNSS) position time-series, and cloud computing (CC) resources. We first briefly describe the basic rationale of the adopted P-SBAS processing approach, tailored to deal with S-1 IWS SAR data and to be implemented in a CC environment, highlighting the innovative solutions that have been introduced in the processing chain we present. They mainly consist in a series of procedures that properly exploit the available GNSS time series with the aim of identifying and filtering out possible residual atmospheric artifacts that may affect the DInSAR measurements. Moreover, significant efforts have been carried out to improve the P-SBAS processing pipeline automation and robustness, which represent crucial issues for interferometric continental scale analysis. Then, a massive experimental analysis is presented. In this case, we exploit: (i) the whole archive of S-1 IWS SAR images acquired over a large portion of Europe, from descending orbits, (ii) the continuous GNSS position time series provided by the Nevada Geodetic Laboratory at the University of Nevada, Reno, USA (UNR-NGL) available for the investigated area, and (iii) the ONDA platform, one of the Copernicus Data and Information Access Services (DIAS). The achieved results demonstrate the capability of the proposed solution to successfully retrieve the DInSAR time series relevant to such a huge area, opening new scenarios for the analysis and interpretation of these ground deformation measurements.

Keywords:

Sentinel-1; DInSAR; P-SBAS; deformation time series; GNSS; DIAS

Graphical Abstract

1. Introduction

The Sentinel-1 (S-1) constellation of the Copernicus Program represents a big revolution within the Earth Observation (EO) scenario, providing an unprecedented operational capability for intensive radar mapping of the Earth surface thanks to its two satellites (Sentinel-1A and B, launched on April 2014 and April 2016, respectively) sharing the same polar-orbit plane and performing C-band synthetic aperture radar (SAR) imaging [1]. In particular, with respect to previous C-band space-borne SAR systems (like those on board the ERS-1/2, ENVISAT, and RADARSAT-1/2 missions [2,3]), the S-1 constellation is characterized by enhanced revisit frequency, spatial coverage, and reliability for operational services and applications requiring long SAR data time series. Moreover, the S-1 constellation archive is available with a free and open access policy, thus easing the data access and enlarging the scientific community interested in its use. Another key element of the S-1 constellation is its main acquisition mode on land, referred to as interferometric wide swath (IWS), which implements the terrain observation by progressive scans (TOPS) technique [4] that guarantees a very large spatial coverage, with a nominal footprint extending for about 250 km, and is specifically designed for interferometric applications. In addition, the S-1 constellation has stringent requirements based on a high attitude and orbit accuracy, and it is intrinsically characterized by small spatial and temporal baselines, with an “orbital tube” of about 200 m nominal diameter and a revisit time of six days (12 days in the case of only one operating satellite).

These characteristics make the S-1 IWS data particularly suitable for exploitation through (conventional and advanced) differential SAR interferometry (DInSAR) techniques [5,6,7], thus opening the possibility to investigate Earth surface deformation phenomena at unprecedented spatial and/or temporal scales. Among the advanced DInSAR techniques, a well-known and widely used one is the small baseline subset (SBAS) approach [8] that investigates surface displacements through the generation of deformation time series and the corresponding mean velocity maps. The SBAS technique has extensively proven its effectiveness to perform analyses with centimeter- to millimeter-level accuracy in several scenarios [2,9,10,11] and is capable of carrying out investigations at different spatial scales and with multi-sensor data [12,13,14,15]. Moreover, a parallel algorithmic solution for the SBAS approach, referred to as the parallel small baseline subset (P-SBAS) technique, was recently developed [16]. It was first applied to process ENVISAT [16,17,18] and COSMO-SkyMed SAR data archives [18] and was subsequently extended to S-1 IWS data [19,20]. In particular, this approach permits the generation of DInSAR time series by taking full benefit of the structure of the S-1 IWS data, which are composed of bursts that can be considered as separate acquisitions and then independently processed up to DInSAR interferogram generation [19,20]. Accordingly, the processing is intrinsically parallelizable with respect to such independent input data portions and, therefore, the P-SBAS approach exploits this coarse granularity parallelization strategy in the majority of the steps of the processing chain [19,20]. Moreover, more sophisticated parallelization approaches are implemented for the steps that are particularly intensive from the computational viewpoint, by exploiting both multi-node and multi-core programming techniques on parallel computing architectures such as cluster, grid, and above all, cloud computing infrastructures [18,21,22].

In this work we present an automatic processing pipeline that exploits cloud computing (CC) resources, based on an advanced version of the P-SBAS processing chain for continental scale deformation time series generation, which benefits from the joint exploitation of S-1 IWS SAR data and global navigation satellite system (GNSS) position time series. Starting from the original P-SBAS processing chain, we describe the implemented innovative solutions that have been introduced. They mainly consist of a series of procedures that benefit from the available GNSS time series, which are properly selected and filtered in order to identify the stations whose signals are not altered by steps, are sufficiently extended in time, and are not affected by local displacement effects, thus accounting for regional (often tectonic) phenomena only. Based on these data, the implemented procedure allows us to identify and filter out possible residual atmospheric artifacts that may affect the DInSAR measurements [23,24]. Moreover, significant efforts have been carried out to improve the P-SBAS pipeline automation and robustness, which represent crucial issues for interferometric continental scale analysis. Subsequently, we present a massive experimental analysis based on the implemented solution. In this case we consider the whole archive of S-1 IWS SAR images and continuous GNSS position time series, the latter provided by the Nevada Geodetic Laboratory at the University of Nevada, Reno, USA (UNR-NGL) [25], collected over a large portion of Europe. Concerning the exploited CC resources, we use the ONDA DIAS platform [26], one of the Copernicus Data and Information Access Services (DIAS), selected through a public tender, for the download and processing of the S-1 images, and for subsequently storing the generated interferometric products. Finally, we present a discussion on the effectiveness of the proposed processing pipeline to successfully retrieve DInSAR time series at continental scale in an automatic, efficient and robust way.

The paper is organized as follows. Section 2 describes the S-1 IWS P-SBAS processing chain and the relevant pipeline implemented within CC environments, highlighting the main automation aspects and algorithmic advancements; the latter have been introduced to effectively exploit the GNSS position time series for the generation of continental scale deformation time series, and are fully detailed in Section 3. Section 4 describes the set of S-1 IWS and GNSS data, as well as the CC resources exploited for the massive experimental interferometric analysis presented, whose results are shown in Section 5. Section 6 is devoted to some conclusive remarks about the developed pipeline and to the new perspectives related to this kind of advanced DInSAR analysis.

2. The S-1 P-SBAS Approach: Basic Rationale and Cloud Computing Implementation

We present in this section a short overview of the exploited S-1 IWS P-SBAS processing chain, whose block diagram is sketched in Figure 1. Note that it represents a modified version of the one reported in [20], since we have introduced an additional step aimed at identifying and filtering out possible residual atmospheric artifacts that may affect the DInSAR measurements. This result is achieved through the proper exploitation of the available GNSS position time series, whose details are provided in the following section. Moreover, we have improved the automation and robustness of the S-1 P-SBAS processing chain from the algorithmic point of view by introducing several check mechanisms to guarantee its fully unsupervised execution.

In the following, we briefly describe the main blocks characterizing this approach and the basic rationale of its implementation in the exploited CC environment, without detailing the implemented algorithmic solutions and the parallelization techniques exploited for each processing step, which are extensively discussed in [19,20]. In this sense, it is worth noting that the black and blue colors used to outline the different blocks shown in Figure 1 allow us to easily discriminate the processing steps that are intrinsically sequential (black) from those where the concurrent jobs are executed in parallel (blue) by distributing the input data processing among different computing nodes. Among the latter, the steps that benefit from the aforementioned burst partitioning granularity, which represents a key element of the parallelization strategy of the exploited S-1 processing chain [19,20], are marked by red color texts.

In the first block of the diagram in Figure 1 (block A) we implement the operations for handling and ingesting the input data, represented by the sequence of S-1 single look complex (SLC) burst images, the orbital information associated to each SAR acquisition, the digital elevation model (DEM) of the investigated zone, and the GNSS position time series available for the area that, in our case, are provided by the Nevada Geodetic Laboratory at the University of Nevada, Reno, USA (UNR-NGL) [25].

Once the input data are correctly ingested, block B performs a series of sequential operations, comprising the selection of the reference SAR geometry for the registration of the entire SAR dataset, the identification of the area of interest, as well as of the burst images to be processed. Moreover, a first list of interferometric data pairs involving the available SAR acquisitions is created, which is then exploited within the subsequent interferograms generation and noise filtering operations [27].

Block C performs the conversion of the available DEM into the radar (azimuth and range) coordinates through the range Doppler equations [7], followed by the estimation, for each pixel of each SAR scene, of the range and azimuth distances with respect to the reference orbital position.

Block D concerns the co-registration of each burst image with respect to the reference acquisition, the burst interferogram generation, and subsequent correction to compensate the effects of possible residual azimuthal misregistration errors. Note also that a spatial multi-look operation is carried out to mitigate the decorrelation noise effects within the generated burst interferograms and to reduce the amount of data to process.

At this stage, in block E the corrected multi-look interferograms (and the corresponding spatial coherence maps) of adjacent bursts are accurately assembled through an interferogram mosaicking operation to generate the differential interferograms of the whole investigated area, which is hereafter referred to as “frame”. Moreover, on these interferograms, the noise filtering procedure discussed in [27] is applied.

Subsequently, a selection, among all the considered small baseline interferometric pairs, of an optimized triangular network within the perpendicular baseline/time plane is computed in block F [27]. In particular, regarding the temporal and spatial baselines, we select a maximum value of 150 days for the former and of 200 m for the latter. This implies that essentially no constraint is applied to the spatial baselines because the orbital tube diameter is about 200 m.

The sequence of the so-identified multi-look noise-filtered small baseline interferograms of the frame is then unwrapped through the extended minimum cost flow (EMCF) phase unwrapping (PhU) algorithm [28] (block G). Then, in block H, a first estimate of the deformation time series is retrieved by applying the SVD method, following the lines of the original SBAS approach [8]. This block also implements a procedure for the compensation of possible topographic phase residuals. Moreover, it includes the estimation and removal of atmospheric artifacts by taking into account that they are typically correlated in space and poorly in time [8], as well as by benefiting from their correlation with topography [29].

Block I represents an additional step implemented within the presented P-SBAS processing chain, which is discussed in details in the following section. In this block, the available GNSS measurements are used to identify and filter out possible residual atmospheric artifacts that may affect the DInSAR measurements [23,24].

Lastly, block J computes the final displacement time series and the corresponding mean deformation velocity map of the investigated frame, generated in a geographic/cartographic reference system.

Concerning the CC implementation of the S-1 IWS P-SBAS processing chain shown in Figure 1, and the corresponding job scheduling, we remark that they are carried out through a fully automatic pipeline able to exploit the available resources, which follows the lines of the solution presented in [19,20] and is sketched in Figure 2. In particular, it starts by downloading and unpacking the S-1 input data in parallel, from the input data archive, by ingesting as inputs the list of S-1 IWS SLC burst images relevant to the investigated frame, and identified by the vertexes of the area of interest to be analyzed. Then, it distributes the S-1 IWS data among the available computing nodes and manages their parallel processing by using both multi-node and multi-core programming techniques. We also stress that different parallelization strategies are adopted throughout several steps of the processing chain of Figure 1, depending on the implemented algorithm structure, as well as on the amount and the type of data they process and on the available computing resources required (RAM, I/O bandwidth, CPUs consumption, etc.) [18,21,22].

At the end of each processing, the developed pipeline performs the upload of the generated DInSAR results on the storage archive dedicated to their preservation. Moreover, as soon as the processing of one frame is completed on a computing node, it launches a new query to the archive in order to proceed with the next S-1 frame. The final results can be visualized and analyzed through appropriate GIS and WebGIS platforms [30].

The overall pipeline is written through a set of bash scripts, which allow exploiting any kind of Linux-based computing environment [19].

Moreover, we remark that in the context of massive DInSAR analyses an effective automation of the input data ingestion and efficient processing, and of the storing of the achieved results represents a crucial point to be addressed. Accordingly, in order to guarantee the fully unsupervised execution of the S-1 P-SBAS processing, we have introduced several check mechanisms at the algorithmic level to guarantee the robustness and the reliability of the pipeline.

In particular, once selected the area of interest relevant to the frame to be processed, there is a first algorithm check verifying that the selected burst images fully cover the selected area, thus avoiding gaps in the DInSAR coverage. Another check mechanism is implemented to investigate if an S-1 input data has been processed more than once with different IPF version [31] in order to select the more recent one. Additional checks are focused on automatically identifying the reference coherent pixel on each frame by also exploiting the available GNSS information or, if unavailable, the triangular coherence [20].

3. Joint Exploitation of GNSS and DInSAR Measurements

The GNSS measurements can provide three-dimensional information, with millimetre accuracy, on the ongoing deformations (East, North, and vertical components) through the accurate evaluation of the receiver stations position. However, the spatial distances among the GNSS stations are often too large to capture localized displacement patterns, while they are typically appropriate for the study of the low spatial frequency deformation fields, often due to regional (tectonic, in several cases) phenomena. On the other hand, the DInSAR techniques provide, for SAR data acquired from a single illumination geometry, one-dimensional surface deformation measurement along the radar Line of Sight (LOS), but with a wide spatial density; this peculiarity makes DInSAR an incomparable technology for reconstructing high spatial frequency deformation patterns. Accordingly, the joint use of GNSS and DInSAR technologies appears to be particularly effective for the retrieval of the deformation signals affecting a specific area of interest [23,24], and is even more important for continental scale analysis like the one presented in this work. Accordingly, this is the key task of the approach discussed in the following.

To clarify the basic rationale of the presented solution, let us first focus on the LOS DInSAR displacements, whose expression can be written as follows:

d_{L o S} = d_{H F} + d_{L F} + η_{A P S} + η_{H F} + η_{L F}

(1)

where

d_{H F}

is the high spatial frequency deformation component due to localized deformation patterns (such as landslides, local subsidence, etc.),

d_{L F}

represents the low spatial frequency deformation relevant to possible regional displacements,

η_{A P S}

is due to the tropospheric noise affecting the DInSAR signals,

η_{H F}

accounts for high spatial frequency deformation effects due to PhU errors, i.e.,

2 π

multiples that were wrongly estimated within the unwrapping operation of the interferometric phase [27,28], and

η_{L F}

is relevant to low spatial frequency artefacts due to atmospheric residuals. It is worth noting that possible DEM errors have been estimated and compensated for in block H although, because the S-1 system is characterized by a very narrow orbital tube [20], the corresponding phase errors can be considered negligible and therefore they are not included in Equation (1). Starting from the

d_{L o S}

measurements, we want to correctly retrieve the

d_{H F}

and

d_{L F}

signals in Equation (1).

Concerning the P-SBAS processing chain, we remark that it already includes methods for the PhU errors mitigation (factor

η_{H F}

, in Equation (1)) [32] and for the APS removal (factor

η_{A P S}

, in Equation (1)) [8]. Accordingly, our main efforts have been focused on the mitigation of the impact of the factor

η_{L F}

(see Equation (1)). In particular, the overall analysis is focused on the identification and removal of the mean rate of the signal component

η_{L F}

, which is referred hereinafter as

v_{η_{L F}}

, through the profitable exploitation of the low pass frequency deformation signals retrieved from the GNSS measurements. To do this, we have developed a rather simple solution, which is depicted in the block diagram of Figure 3 and discussed in detail as follows.

3.1. GNSS Data Screening and Selection

The starting point of our algorithm is a preliminary screening of the measurements accessible through the available GNSS stations. In particular, this allows the identification of the GNSS stations whose measurements can be effectively used to filter out the signal component

v_{η_{L F}}

To do this, as detailed in the following, it is essential to identify the stations whose signals are:

Sufficiently extended in time (in our case at least 2 years) and “cleaned” of possible jumps and artefacts, which may affect the correct estimation of the mean deformation rate components;
Affected by regional scale signals only, which are therefore spatially correlated and that do not account for high spatial frequency deformation signals (i.e., localized displacements).

To clarify the need of the above-mentioned data screening, we remark that the GNSS stations time series can sometimes be perturbed by the presence of gaps in the data acquisitions and by steps, which can be generated by different causes, both anthropogenic and natural. In particular, software modifications or physical changes of the antenna receiver, as well as the occurrence of a seismic event, are the main reasons for the presence of these steps within the time series that, unfortunately, are not always identified through the auxiliary files, where the dates in correspondence of these issues are highlighted.

Accordingly, our data screening step identifies first the stations whose signals have a significant temporal extension that, in our case, corresponds to at least two years. Then, from each of them we extract the time interval that, with respect to the presence of possible steps in the time series, is reliable concerning the evaluation of the mean rate component. Accordingly, the GNSS stations whose signals are not consistent with the above characteristics are rejected from the procedure.

In Figure 4, five examples of the GNSS stations time series selection are shown. In particular, the plots on the left side show the original time series (we consider the East/West component as provided by the UNR-NGL web-site (http://geodesy.unr.edu)), while on the right, the corresponding ones following the automatic data screening procedure, are represented. Note that the plots 1 and 2 are representative of GNSS stations affected by seismic events: in particular, these two stations show the steps due to the 2016 Central Italy seismic sequence (24 August 2016, 30 October 2016, and 18 January 2017). Plots 3 and 4 outline some anthropic operations performed on the GNSS stations, which in some cases corrupts the GNSS station deformation rate estimation. Finally, we show that plot 5, although not supported by the auxiliary information, is filtered by the algorithm, which identifies the large temporal gap in the data and rejects the shorter part of the available time series, thus allowing a proper mean deformation rate estimation.

At this stage, the implemented procedure identifies the GNSS stations affected by regional signals only. To do this we exclude the GNSS stations characterized by a significant vertical deformation velocity (greater than 0.5 cm/year), because this is not consistent with long-term regional deformation signals and we take advantage of the strong spatial correlation of these phenomena, whose planar component velocities are spatially correlated for tens of squared kilometres [33]. In particular, we evaluate a similarity rate index of the GNSS stations that fall within a specific spatial window that in our case is empirically set to ~30 × 30 km². This similarity index, which we estimate for each GNSS station, is expressed as follows:

R_{i} = \frac{1}{n - 1} \sum_{j = 1}^{n - 1} [\frac{v_{i} v_{j}}{v_{i}^{2}}] = \frac{1}{n - 1} \sum_{j = 1}^{n - 1} [\frac{v_{j}}{v_{i}}]

(2)

where n represents the total number of the GNSS stations that fall within the above-mentioned spatial window,

v

is the GNSS mean displacement rate, while i is the station where we are calculating the similarity index, and j is the index relevant to the remaining n × 1 stations that we compare with i. The spatial similarity rate index is designed to be equal to 1 when the deformation velocities of the stations falling within the spatial window are the same, while the similarity becomes worse for values of this index that move away from 1.

The evaluation of the similarity index in Equation (2) is typically carried out only on the planar components of the GNSS-based estimate of the mean deformation rate, which are expressed as follows for the ith GNSS station:

R_{e w_{i}} = \frac{1}{n - 1} \sum_{j = 1}^{n - 1} [\frac{v_{e w_{i}} v_{e w_{j}}}{v_{e w_{i}}^{2}}] = \frac{1}{n - 1} \sum_{j = 1}^{n - 1} [\frac{v_{e w_{j}}}{v_{e w_{i}}}]

(3)

R_{n s_{i}} = \frac{1}{n - 1} \sum_{j = 1}^{n - 1} [\frac{v_{n s_{i}} v_{n s_{j}}}{v_{n s_{i}}^{2}}] = \frac{1}{n - 1} \sum_{j = 1}^{n - 1} [\frac{v_{n s_{j}}}{v_{n s_{i}}}]

(4)

In particular, only the GNSS stations for which

(R_{e w_{i}}, R_{n s_{i}}) \in [0.95, 1.05]

are considered, within the DInSAR deformation signals refinement step discussed in the following paragraph.

In Figure 5, we represent the planar components (East-West and North-South) of the GNSS stations of the NGL catalogue, located in the entire European territory under investigation in this work. It is evident from the visual inspection of Figure 5 that the stations selected for our processing are highly spatially correlated and with no evident singularities present.

3.2. DInSAR Deformation Signals Refinement

In this paragraph, we describe the last steps of our procedure for the refinement of the DInSAR deformation signals relevant to each investigated SAR frame. In particular, this procedure aims to identify and remove possible low spatial frequency displacements artefacts, mainly due to residual atmospheric phase errors.

In Figure 6, we show a panel describing, step by step, the implemented procedure applied to the P-SBAS results generated for a selected frame that is relevant to a large portion of the Central-Northern Italy regions, included in the processing experiment presented in Section 5. Note that we select this frame because it is affected by different displacement phenomena, thus highlighting the efficacy of the implemented algorithm. In particular, the included Po Valley is historically affected by intense subsidence phenomena, as well as the city of Bologna [34]. Moreover, in the same frame we can appreciate the well-known displacements relevant to the Firenze–Prato–Pistoia basin [35,36], the subsidence phenomena affecting the Northern Adriatic Sea coastline [37], and in the South-East corner of the frame the impact of the August–October 2016 seismic sequence [38,39].

Figure 6A shows the starting point for the developed procedure. In particular it represents the mean displacement velocity map obtained from the original P-SBAS approach (see block H in Figure 1). Figure 6B depicts the low spatial frequency deformation pattern retrieved from the interpolation (Kriging interpolation [40]) of the GNSS stations planar displacement velocities projected along the satellite radar LOS. Clearly, for this step, we exploit only the GNSS stations returned by the screening operation and for values of the similarity rate index close to 1, as described in Section 3.1.

The mean deformation velocity retrieved by subtracting the GNSS-derived regional deformation pattern (Figure 6B) from the starting velocity (Figure 6A) corresponds to the input for the velocity calibration procedure. Such a procedure is equivalent to apply a simple low pass spatial filter for the estimation and removal of possible artefacts. Figure 6C shows the result of the calibration step that allows us to correct the retrieved deformation signals.

At this stage, to reconstruct the overall displacement field affecting the area of interest, we can add the low and high spatial frequency deformation components (Figure 6B,C, respectively). The mean deformation velocity map corresponding to the final result is shown in Figure 6D. Moreover, the above-described procedure is also graphically summarized with a block diagram at the bottom of Figure 6.

As an additional remark, we underscore that both results shown in Figure 6C,D can be separately exploited for different scientific purposes, thus increasing the range of the value-added results offered by these advanced DInSAR analyses.

We finally present the results of the comparison between the refined DInSAR deformation time series and those of the available GNSS stations relevant to the frame of Figure 6. Obviously, we exploit in this comparison only the GNSS stations that have not been accounted within the DInSAR signals refinement step.

In Figure 7 we show the DEM of the area, highlighting the sites of all the GNSS stations that have been investigated in our study; they are identified by blue and yellow triangles, representing the stations exploited for and excluded from the DInSAR deformation signals refinement step, respectively. In Figure 8, we report a selection of plots showing the comparisons between the deformation time series retrieved through the developed S-1 IWS P-SBAS processing pipeline (black triangles) and the corresponding LOS-projected GNSS ones (red stars) for the stations labelled in Figure 7 with capital letters from A to R. We remark that, since the GNSS time series themselves are affected by noise (although they are assumed as reference), we considered a smoothed version of the GNSS time series in our comparative analysis.

Figure 8A–O show the plots relevant to some GNSS stations that were excluded from the DInSAR deformation signals refinement procedure. The good agreement between the P-SBAS and LOS-projected GNSS measurements is evident, as also testified by the standard deviation values of the difference between the two time series, computed in the temporal window common to both measurements, ranging from 0.23 and 0.49 cm. Moreover, we also show in Figure 8P–R the plots for three GNSS stations where, for different reasons, the performed comparison is not satisfactory. In particular, Figure 8P is relevant to a GNSS station located in the area that was strongly struck by the 2016 Central Italy seismic sequence. In this case, the discontinuity characterizing the DInSAR time series in correspondence to the earthquake is clearly underestimated due to PhU errors [41,42]. In addition, Figure 8Q,R present the plots of two GNSS stations where the effects of the antenna code changes, as reported in the GNSS auxiliary files, are evident.

4. Exploited SAR and GNSS Data and Cloud Computing Resources

We describe in this section the input SAR and GNSS data we exploit in this study. Moreover, we briefly describe also the cloud computing resources and architecture we use to perform the whole P-SBAS processing.

4.1. SAR and GNSS Data

As stated already, the SAR images we analyze have been acquired by the Sentinel-1A and Sentinel-1B satellites. The overall investigated area has been divided in 175 frames, from 21 descending orbit tracks, with an average size of almost 200 × 250 km (azimuth and range, respectively). An overlap between frames has been preserved, covering almost 20% of the frame size for two subsequent ones in the azimuth direction. This overlapping area allows us to perform consistency tests on the analyzed zone, in order to guarantee an additional check for the continuity preservation of the achieved result. In Figure 9, we show the footprint of the 175 S-1 frames processed in this study.

The time series of every frame have been retrieved by using data acquired from March 2015 to September 2018 with the average number of SAR images for each frame varying between 400 and 450, depending on the different regions. Moreover, for the frames relevant to Italy, our analysis has been extended to March 2020, since this study also supports a project led by the Italian Ministry of Economic Development aimed at investigating the deformation phenomena affecting the Italian territory.

Globally we have used almost 72,000 SAR images covering a whole area of about 4,500,000 km², including the overlap among the considered frames.

Concerning the GNSS stations, we have used, as previously mentioned, the data provided by UNR-NGL [25] in the framework of the “Plug and Play GPS” project that takes the raw GNSS data from more than 17,000 stations all around the world and makes publicly available via a web-site (http://geodesy.unr.edu) the data products (including metadata, plots of position coordinates, etc.) for all the geodetic quality GNSS stations. The starting number of GNSS stations, whose data have been downloaded and analyzed in this study, is 4227. Following the data screening and selection procedure described in Section 3, we have exploited 2477 stations relevant to the investigated area. In Figure 10, we show the position of the stations: blue markers represent the stations used after the screening, while the yellow ones represent the stations whose data have been neglected because they are considered unsuitable for the implemented DInSAR deformation signals refinement operation.

4.2. Cloud Computing Resources

As anticipated before, we perform the CC implementation of the S-1 IWS P-SBAS processing pipeline discussed in Section 2, by exploiting the computing resources provided by the ONDA DIAS CC platform [26]. In particular, we use six computing nodes, each one equipped according to the characteristics reported in Table 1. Note that the six nodes are connected to the S-1 IWS data archive, which is located on an object storage, via NFS with a guaranteed download bandwidth of 2 Gb/sec (see Figure 11). As mentioned, having the data archive close to the computing resources and accessible through a dedicated and fast connection is a key requirement when the processing of such a huge data volume is envisaged, because of the data transfer overhead. In this framework we also remark that we fully benefit from the specific APIs made available by ONDA in order to effectively and automatically perform the S-1 images query and download operations.

5. Experimental Results

We processed through the CC-based S-1 P-SBAS pipeline described in Section 2 the S-1 IWS data-set acquired over the wide portion of the European territory detailed in Section 4.1. As stated before, the data-set considered is composed of 175 frames (see Figure 9) and for each of these, the selected SAR image temporal sequences were individually processed. We remark that no spatial baseline constraint was imposed in the interferometric pair selection exploited in our processing thanks to the narrow orbital tube characterizing the S-1 constellation [20]. Moreover, we exploited the 1-arcsec SRTM DEM to generate the analyzed DInSAR interferograms on which we performed a complex multi-look operation with twenty looks in the range direction and five looks in the azimuth one, thus obtaining a pixel dimension of about 80 × 80 m. The same DEM was also used to geocode the computed deformation time series and the corresponding mean velocity maps.

In Figure 12 we present the mosaicking of the obtained mean deformation velocity maps of the overall investigated area, geocoded and superimposed on an optical image. The displayed overall mean velocity map provides information on the ongoing deformation mean values of the coherent pixels, identified by considering those with a temporal coherence value [28] greater than a selected threshold that in our analysis is set equal to 0.9 for all the analyzed data-sets. Overall, we detect ~120,000,000 coherent pixels for which the corresponding deformation time series are available. Moreover, the zones exhibiting lower temporal coherence values are properly masked out. The velocity maps in Figure 12 show a large variety of localized deformation patterns due to different phenomena (volcanic activities, earthquakes, subsidence phenomena, slope instabilities, man-made activities, etc.), with the most evident one being the well-known deformation pattern associated with the seismic crisis that struck Central Italy during the August–October 2016 period, already extensively discussed in [38,39].

Starting from the results of Figure 12, and to provide a clearer idea about the generated products, we present in the following the zoom-in views of the mean deformation velocity map relevant to six areas characterized by significant ground displacements related to natural phenomena (Figure 13, Figure 14 and Figure 15) and man-induced activities (Figure 16, Figure 17 and Figure 18), which are outlined by white rectangles. Moreover, for each area, we show the plots of the deformation time series of some selected coherent pixels.

In particular, in Figure 13a, we show the zoom-in view of the mean deformation velocity map portion relevant to the Campi Flegrei caldera (Southern Italy), an active volcanic area located West of the city of Napoli. It is characterized by the bradyseism phenomenon consisting in a slow lifting and lowering movement of the ground; since late 2005 [43,44,45,46,47] the caldera has been experiencing an almost continuous uplift of the ground surface. Indeed, the map in Figure 13a shows a spatially extended deformation pattern involving the whole caldera with a maximum LOS deformation rate of about 7.5 cm/year. We also present the plot (Figure 13b) of the displacement time series of a pixel located in the maximum deformation area, identified by a white star and labeled as P1 in Figure 13a, characterized by a cumulative displacement of about 37 cm. The temporal evolution of the displacement for the selected pixel shows the deformation behavior characterizing the caldera since 2015. As a matter of fact, we observe a rather continuous uplift typified by different deformation rates within the whole period of observation. In addition, a decrease of the deformation rate, approximately spanning from mid-2016 to mid-2017, is also clearly visible. Finally, in Figure 13c, we report the mean deformation velocity values for the coherent pixels deployed along the section crossing the maximum uplifting zone and identified in Figure 13a by a dashed white line. It is evident that the caldera shows a general uplift characterized by an almost radial symmetry with respect to the maximum displacement located in the Rione Terra area (P1 in Figure 13a).

In Figure 14a, we highlight a sketch of the mean deformation velocity map relevant to the Lefkada Island (Greece). It outlines a spatially extended deformation pattern, characterized by both negative and positive LOS-projected displacement signals, which are related to the Mw 6.5 earthquake that struck the island the 17 November 2015 [48,49]. Moreover, we present the corresponding displacement time series relevant to two selected pixels (labeled as P1 and P2 in Figure 14a and marked by white stars) located in the maximum co-seismic deformation zone. They clearly show the deformation signal associated with the seismic event, see the discontinuity characterizing the time series of Figure 14b,c in correspondence to the earthquake (highlighted by a vertical dashed red line).

In Figure 15, we show a zoomed-in view of the mean deformation velocity map relevant to the Northern Black Sea coast, between the Varna and Kavarna towns (Bulgaria). This area is characterized by numerous landslides, mainly due to soil over-moistening, marine erosion, slope cutting and seismic impacts, which have caused serious damage to buildings and infrastructure over time [50]. Figure 15a reports the LOS mean deformation velocity map of the area while Figure 15b,c display the corresponding displacement time series relevant to two selected pixels (labeled as P1 and P2 in Figure 15a and marked by white stars), located in Balchik town and in Kranevo village (both between Varna city and Kavarna town), respectively. These plots clearly show the ground movements, characterized by different rates (about 1 and 1.8 cm/year for P1 and P2, respectively), related to the landslides.

Moving from the ground deformation analysis related to natural phenomena to human-induced activities, we benefit first from the DInSAR capability to investigate displacements caused by mining activities [51,52]. In particular, we show in Figure 16a,b two zooms of the map of Figure 12 relevant to rather wide zones located in Germany (a) and Poland (b), respectively, where significant mining activities are carried out. In Figure 16c, we report a zoomed-in view in correspondence to the Nochten coalmine and surroundings, located in the Saxony-Anhalt State (Germany), where the effects of the sub-surface exploitation activities on the ground deformation are evident. In addition, we show in Figure 16e the displacement time series relevant to a selected pixel (labeled as P1 in Figure 16c and marked by a white star) that reveals cumulative displacements of more than 45 cm in the considered time period. Similarly, in Figure 16d we show a zoomed-in view of the Katowice coalmines, located in the Silesian Region (Southern Poland). As for the German case, the ground deformation related to mining actions are clear. The plot of the displacement time series (see Figure 16f) relevant to a selected pixel (labeled as P2 in Figure 16d and marked by a white star) located in the area shows a clear non-linear deformation behavior, which is effectively retrieved through our P-SBAS-based analysis.

We remark that, following the deformation time series retrieval operation, the S-1 IWS P-SBAS processing chain automatically supplies various value-added information on the ongoing deformation phenomena. Indeed, in addition to the estimate of the mean deformation velocity, it also automatically identifies possible deformation signals characterized by a periodic behavior with respect to time, such as those relevant to oil and gas extraction/storage activities and water withdrawal, which represent significant DInSAR applications in the energy source storage [53] and hydrogeological [54] fields. This is done by computing the correlation coefficient of each pixel’s detrended time series with a sinusoid with a one-year period. In Figure 17a we show the peak-to-peak amplitude of the sinusoid signal in an area located in Northern Italy, close to the city of Milano. Note that pixels with peak-to-peak amplitude values equal to 0 are those for which the correlation coefficient is smaller than 0.8. The oscillation amplitude map clearly identifies three distinct areas that correspond to three gas storage sites located in the zone. In Figure 17b,c we show the displacement time series relevant to three selected pixels, labeled as P1, P2, and P3 in Figure 17a and marked by black stars, located in correspondence to the Settala, Sergnano, and Ripalta Guerina sites, respectively. The displayed time series reveal periodic deformation signals that correspond to the human-induced actions carried out to ensure gas supply during the winter. Indeed, the plots clearly show the effect of the injection of gas into storage during periods of low demand (approximately from March to October), which produces a raising of the soil over the reservoir area. Conversely, it is evident that a subsidence effect during periods of peak demand (approximately from November to February) occurs due to the withdrawal of gas migration from storage into the town’s gas supply system. As a further example of this kind of analysis, we show in Figure 18 the oscillation amplitude map relevant to the Firenze-Prato-Pistoia basin (North Italy). In this case, the displacement time series relevant to a selected pixel (labeled as P1 in Figure 18a and marked by a black star) shows a clear subsidence signal, caused by the large amount of ground-water extracted due to the industrial activities carried out in this area [35,36]. Moreover, this subsidence trend is modulated by a seasonal oscillation, due to the groundwater withdrawal and recharge during the summer and winter seasons respectively, see Figure 18b.

It is also worth noting that the two examples reported in Figure 17 and Figure 18 identify deformation signals characterized by oscillations that peak in different periods: approximately in November for the first case and February for the second one. This temporal peak location is a key parameter to discriminate the occurrence of such different displacement phenomena.

6. Conclusions

The amount and availability of remote sensing satellite data have exponentially grown in recent years, significantly contributing to the on-going big data scenario. In particular, a key role is played by the Copernicus program of the European Union through the Sentinel missions, characterized by its free and open data policy [55]. In this framework, the wide spatial footprint, the high temporal acquisition rate, the accurate synchronization among subsequent acquisitions, the selection of a unique acquisition mode over land, and the accurate orbital information with a small tube diameter, which characterizes the Sentinel-1 constellation [1], have already represented a revolution for the possibility to operationally apply the DInSAR techniques to investigate Earth surface deformation phenomena at unprecedented spatial and/or temporal scales. In this paper, we present a CC-based solution to perform, through the P-SBAS processing chain, advanced DInSAR analyses from S-1 IWS SAR images and GNSS position time series at continental scale in an automatic, efficient, and robust way.

To achieve such a task, processing and storing huge data flows, including GNSS, DEM, and orbital information, in addition to very large S-1 IWS SAR data volumes, are required. To do this, there is a need, on the one hand, for proper DInSAR processing algorithms and, on the other hand, for appropriate computing architectures, wide storage capacity, and high performance processing resources. In this framework, key problems to be addressed are also those relevant to the automation level and the reliability of the entire processing strategy. Indeed, carrying out continental scale DInSAR analyses results in the ability to process large volumes of SAR data, preventing algorithmic errors, processing and computing failures, and hardware faults as much as possible.

From the algorithmic point of view the developed S-1 IWS P-SBAS processing chain has been enriched with a processing step that, starting from the available GNSS measurements, following a proper screening and selection operation, allows identifying and filtering out possible residual atmospheric artifacts that may affect the DInSAR measurements, thus generating refined DInSAR products (see Section 3). Moreover, regarding the automation and the robustness of the processing chain and of the relevant CC-based pipeline, several check mechanisms have been introduced in order to guarantee a fully unsupervised processing execution, as thoroughly discussed in Section 2.

With respect to the architectural issues, the main aspects impacting the efficiency of the S-1 P-SBAS pipeline are: i) proximity of data to the HPC infrastructure, ii) large storage capacity, and iii) wide computing capability in order to handle the requested massive parallel computation. Therefore, cloud computing, whose main characteristics are high availability, reliability, and robustness, has revealed as a natural choice to implement the presented continental scale S-1 IWS P-SBAS processing pipeline. In this case, we have successfully demonstrated that the European DIAS Cloud Services, in particular the exploited ONDA Cloud Platform, are suitable to host our fully automated pipeline for downloading and processing the analyzed S-1 images, and for storing the final S-1 interferometric products.

We further remark that, through our processing experiment, we have also successfully verified the high reliability and compliance with respect to the defined formats of the S-1 images archives, which represents another very relevant element supporting the possibility to operationally carrying out such massive DInSAR data processing.

The presented results, obtained through an extensive experimental analysis on a large portion of Europe (see Section 5), clearly show that the implemented CC-based S-1 P-SBAS pipeline has proved to be successfully suitable for massive and automatic continental scale DInSAR processing, which has been completed in six months and at relatively low computing costs amounting to about € 60k. In particular, the presented pipeline allowed the production, for the presented experiment, of the deformation time series and the corresponding velocity maps relevant to each single frame in six days, on average. Nevertheless, we underline that possible improvements of the proposed solution can already be envisioned. For instance, we report that an enhancement for improving the storage connection and performance could be very valuable. In addition to the above consideration, we remark that although the presented analysis has been carried out by exploiting averaged (multi-look) interferograms, it can be efficiently extended to the full resolution spatial scale by using graphics processing units (GPUs), following the lines of the approach presented in [56]. Moreover, algorithmic improvements concerning the compensation of residual phase unwrapping errors [32] and the integration of atmospheric information [57] for improving the tropospheric artifacts removal can also be very important.

Generating updatable interferometric products at a continental scale provides the possibility to access an unprecedented amount of results containing value-added information that is potentially groundbreaking. For instance, an extensive analysis on the interferometric coherence, performed on a country by country basis for land cover classification and mapping, or vegetation density detection, or the study of the statistics of the retrieved DInSAR deformation time series (for instance, standard deviation or root mean square values) and of their correlation with specific (linear, logarithmic, exponential, periodic, etc.) signals, to identify different deformation phenomena or potential changes that occurred during the observation period, could represent very relevant issues. Moreover, the availability of spatially dense deformation time series over such large areas paves the way to new kind of analyses and interpretation of such data. For instance, the idea of applying machine and deep learning techniques [58,59] to follow the behavior of deformation signals and to automatically extract significant patterns is becoming a very relevant issue but other innovative methodologies could be devised for deeply exploiting these “new generation” sources of DInSAR measurements.

As a final remark, we underline that the results we achieved can be particularly relevant for the Solid Earth community that promoted this experiment as a key study of the EPOSAR service included in the Satellite Data Thematic Core Service of the European Plate Observing System (EPOS) (https://www.epos-ip.org/). Moreover, our findings may also have a positive impact on the future development of the European Ground Motion Service [60], which is expected to start in the near future.

Author Contributions

Conceptualization, R.L., M.B., F.C., C.D.L., M.M. (Michele Manunta), M.M. (Mariarosaria Manzo), G.O., and I.Z.; methodology, R.L., M.B., F.C., C.D.L., M.M. (Michele Manunta), M.M. (Mariarosaria Manzo), G.O., and I.Z.; validation, M.B., C.D.L., and M.M. (Mariarosaria Manzo); data curation, M.B., M.M. (Michele Manunta), and G.O.; writing—original draft preparation, R.L., C.D.L., M.M. (Mariarosaria Manzo), G.O., and I.Z.; visualization, F.C., G.O., and I.Z.; supervision, R.L. and F.C. All authors have read and agreed to the published version of the manuscript.

Funding

This research was funded by the Italian Ministry of Economic Development, the Italian Civil Protection Department, the I-AMICA (PONa3_00363) project, and the EPOS-SP, EOSC-hub, ENVRI-FAIR and OpenAIRE-advance projects of the European Union Horizon 2020 for Research and Innovation Program under Grant 871121, Grant 777536, Grant 824068, and Grant 777541, respectively.

Acknowledgments

This work is supported by the 2018-2019 IREA-CNR/DGS-UNMIG and the 2019-2021 IREA-CNR and Italian Civil Protection Department agreements. On the other hand, we remark that the contents of this paper represent the authors’ ideas and do not necessarily correspond to the official opinion and policies of the Italian Ministry of Economic Development and Italian Civil Protection Department. This study is also supported by the I-AMICA (PONa3_00363) project and the EPOS-SP, EOSC-hub, ENVRI-FAIR and OpenAIRE-advance projects of the European Union Horizon 2020 for Research and Innovation Program under Grant 871121, Grant 777536, Grant 824068, and Grant 777541, respectively. The Sentinel-1 data have been provided through the Copernicus Program of the European Union. The DEM of the investigated zone was acquired through the SRTM archive. The authors thank Geoffrey Blewitt of UNR for his support in the access to the exploited GNSS data, Enrico Priolo of OGS for the interesting discussions on the possible use of the presented results for monitoring oil and gas extraction/storage activities and Paul R. Lundgren of JPL, California Institute of Technology, for proofreading the paper. Finally, the authors thank S. Guarino, F. Parisi and M.C. Rasulo for their support.

Conflicts of Interest

The Authors declare no conflict of interest.

References

Torres, R.; Snoeij, P.; Geudtner, D.; Bibby, D.; Davidson, M.; Attema, E.; Potin, P.; Rommen, B.; Floury, N.; Brown, M.; et al. GMES Sentinel-1 mission. Remote Sens. Environ. 2012, 120, 9–24. [Google Scholar] [CrossRef]
Bonano, M.; Manunta, M.; Pepe, A.; Paglia, L.; Lanari, R. From Previous C-Band to New X-Band SAR Systems: Assessment of the DInSAR Mapping Improvement for Deformation Time-Series Retrieval in Urban Areas. IEEE Trans. Geosci. Remote Sens. 2013, 51, 1973–1984. [Google Scholar] [CrossRef]
Sansosti, E.; Berardino, P.; Bonano, M.; Calo, F.; Castaldo, R.; Casu, F.; Manunta, M.; Manzo, M.; Pepe, A.; Pepe, S.; et al. How second generation SAR systems are impacting the analysis of ground deformation. Int. J. Appl. Earth Obs. Geoinf. 2014, 28, 1–11. [Google Scholar] [CrossRef] [Green Version]
De Zan, F.; Guarnieri, A.M. TOPSAR: Terrain Observation by Progressive Scans. IEEE Trans. Geosci. Remote Sens. 2006, 44, 2352–2360. [Google Scholar] [CrossRef]
Massonnet, D.; Feigl, K.L. Radar interferometry and its application to changes in the Earth’s surface. Rev. Geophys. 1998, 36, 441–500. [Google Scholar] [CrossRef] [Green Version]
Bürgmann, R.; Rosen, P.A.; Fieldingi, E. Synthetic Aperture Radar Interferometry to Measure Earth’s Surface Topography and Its Deformation. Annu. Rev. Earth Planet. Sci. 2000, 28, 169–209. [Google Scholar] [CrossRef]
Franceschetti, G.; Lanari, R. Synthetic Aperture Radar Processing; CRC Press: Boca Raton, FL, USA, 1999. [Google Scholar]
Berardino, P.; Fornaro, G.; Lanari, R.; Sansosti, E. A new algorithm for surface deformation monitoring based on small baseline differential SAR interferograms. IEEE Trans. Geosci. Remote Sens. 2002, 40, 2375–2383. [Google Scholar] [CrossRef] [Green Version]
Casu, F.; Manzo, M.; Lanari, R. A quantitative assessment of the SBAS algorithm performance for surface deformation retrieval from DInSAR data. Remote Sens. Environ. 2006, 102, 195–210. [Google Scholar] [CrossRef]
Lanari, R.; Casu, F.; Manzo, M.; Lundgren, P. Application of the SBAS-DInSAR technique to fault creep: A case study of the Hayward fault, California. Remote Sens. Environ. 2007, 109, 20–28. [Google Scholar] [CrossRef]
Manzo, M.; Fialko, Y.; Casu, F.; Pepe, A.; Lanari, R. A Quantitative Assessment of DInSAR Measurements of Interseismic Deformation: The Southern San Andreas Fault Case Study. Pure Appl. Geophys. 2012, 169, 1463–1482. [Google Scholar] [CrossRef]
Lanari, R.; Mora, O.; Manunta, M.; Mallorqui, J.J.; Berardino, P.; Sansosti, E. A small-baseline approach for investigating deformations on full-resolution differential SAR interferograms. IEEE Trans. Geosci. Remote Sens. 2004, 42, 1377–1386. [Google Scholar] [CrossRef]
Manunta, M.; Marsella, M.; Zeni, G.; Sciotti, M.; Atzori, S.; Lanari, R. Two-scale surface deformation analysis using the SBAS-DInSAR technique: A case study of the city of Rome, Italy. Int. J. Remote Sens. 2008, 29, 1665–1684. [Google Scholar] [CrossRef]
Pepe, A.; Sansosti, E.; Berardino, P.; Lanari, R. On the Generation of ERS/ENVISAT DInSAR Time-Series Via the SBAS Technique. IEEE Geosci. Remote Sens. Lett. 2005, 2, 265–269. [Google Scholar] [CrossRef]
Bonano, M.; Manunta, M.; Marsella, M.; Lanari, R. Long-term ERS/ENVISAT deformation time-series generation at full spatial resolution via the extended SBAS technique. Int. J. Remote Sens. 2012, 33, 4756–4783. [Google Scholar] [CrossRef]
Casu, F.; Elefante, S.; Imperatore, P.; Zinno, I.; Manunta, M.; De Luca, C.; Lanari, R. SBAS-DInSAR Parallel Processing for Deformation Time-Series Computation. IEEE J. Sel. Top. Appl. Earth Obs. Remote Sens. 2014, 7, 3285–3296. [Google Scholar] [CrossRef]
De Luca, C.; Zinno, I.; Manunta, M.; Lanari, R.; Casu, F. Large areas surface deformation analysis through a cloud computing P-SBAS approach for massive processing of DInSAR time series. Remote Sens. Environ. 2017, 202, 3–17. [Google Scholar] [CrossRef]
Zinno, I.; Casu, F.; De Luca, C.; Elefante, S.; Lanari, R.; Manunta, M. A Cloud Computing Solution for the Efficient Implementation of the P-SBAS DInSAR Approach. IEEE J. Sel. Top. Appl. Earth Obs. Remote Sens. 2017, 10, 802–817. [Google Scholar] [CrossRef]
Zinno, I.; Bonano, M.; Buonanno, S.; Casu, F.; De Luca, C.; Manunta, M.; Manzo, M.; Lanari, R. National Scale Surface Deformation Time Series Generation through Advanced DInSAR Processing of Sentinel-1 Data within a Cloud Computing Environment. IEEE Trans. Big Data 2018, 4, 1. [Google Scholar] [CrossRef] [Green Version]
Manunta, M.; De Luca, C.; Zinno, I.; Casu, F.; Manzo, M.; Bonano, M.; Fusco, A.; Pepe, A.; Onorato, G.; Berardino, P.; et al. The Parallel SBAS Approach for Sentinel-1 Interferometric Wide Swath Deformation Time-Series Generation: Algorithm Description and Products Quality Assessment. IEEE Trans. Geosci. Remote Sens. 2019, 57, 6259–6281. [Google Scholar] [CrossRef]
Zinno, I.; Elefante, S.; Mossucca, L.; De Luca, C.; Manunta, M.; Terzo, O.; Lanari, R.; Casu, F. A First Assessment of the P-SBAS DInSAR Algorithm Performances Within a Cloud Computing Environment. IEEE J. Sel. Top. Appl. Earth Obs. Remote Sens. 2015, 8, 4675–4686. [Google Scholar] [CrossRef]
Zinno, I.; Mossucca, L.; Elefante, S.; De Luca, C.; Casola, V.; Terzo, O.; Casu, F.; Lanari, R. Cloud Computing for Earth Surface Deformation Analysis via Spaceborne Radar Imaging: A Case Study. IEEE Trans. Cloud Comput. 2015, 4, 1. [Google Scholar] [CrossRef]
Parizzi, A.; González, F.R.; Brcic, R. A Covariance-Based Approach to Merging InSAR and GNSS Displacement Rate Measurements. Remote Sens. 2020, 12, 300. [Google Scholar] [CrossRef] [Green Version]
Bischoff, C.A.; Ferretti, A.; Novali, F.; Uttini, A.; Giannico, C.; Meloni, F. Nationwide deformation monitoring with SqueeSAR^® using Sentinel-1 data. In Proceedings of the Proceedings of the International Association of Hydrological Sciences; Copernicus GmbH: Göttingen, Germany, 2020; Volome 382, pp. 31–37. [Google Scholar]
Blewitt, G.; Hammond, W.C. Harnessing the GPS Data Explosion for Interdisciplinary Science. Eos 2018, 99, 99. [Google Scholar] [CrossRef]
Nbsp. Available online: https://www.onda-dias.eu/cms/ (accessed on 10 July 2019).
Pepe, A.; Yang, Y.; Manzo, M.; Lanari, R. Improved EMCF-SBAS Processing Chain Based on Advanced Techniques for the Noise-Filtering and Selection of Small Baseline Multi-Look DInSAR Interferograms. IEEE Trans. Geosci. Remote Sens. 2015, 53, 4394–4417. [Google Scholar] [CrossRef]
Pepe, A.; Lanari, R. On the Extension of the Minimum Cost Flow Algorithm for Phase Unwrapping of Multitemporal Differential SAR Interferograms. IEEE Trans. Geosci. Remote Sens. 2006, 44, 2374–2383. [Google Scholar] [CrossRef]
Shirzaei, M.; Bürgmann, R. Topography correlated atmospheric delay correction in radar interferometry using wavelet transforms. Geophys. Res. Lett. 2012, 39. [Google Scholar] [CrossRef]
Buonanno, S.; Zeni, G.; Fusco, A.; Manunta, M.; Marsella, M.; Carrara, P.; Lanari, R. A GeoNode-Based Platform for an Effective Exploitation of Advanced DInSAR Measurements. Remote Sens. 2019, 11, 2133. [Google Scholar] [CrossRef] [Green Version]
The Sentinel-1 IPF Auxiliary Product Specification. Available online: https://sentinel.esa.int/web/sentinel/user-guides/sentinel-1-sar/document-library/-/asset_publisher/1dO7RF5fJMbd/content/sentinel-1-ipf-auxiliary-product-specification;jsessionid=DA0DD0EF437909E50BA09E8F1FD22CD5.jvm1?redirect=https%3A%2F%2Fsentinel.esa.int%2Fweb%2Fsentinel%2Fuser-guides%2Fsentinel-1-sar%2Fdocument-library%3Bjsessionid%3DDA0DD0EF437909E50BA09E8F1FD22CD5.jvm1%3Fp_p_id%3D101_INSTANCE_1dO7RF5fJMbd%26p_p_lifecycle%3D0%26p_p_state%3Dnormal%26p_p_mode%3Dview%26p_p_col_id%3Dcolumn-1%26p_p_col_count%3D1 (accessed on 10 September 2020).
De Luca, C.; Onorato, G.; Casu, F.; Lanari, R.; Manunta, M. A Genetic Algorithm for Phase Unwrapping Errors Correction in the SBAS-DInSAR Approach. In Proceedings of the IGARSS 2019—2019 IEEE International Geoscience and Remote Sensing Symposium; Institute of Electrical and Electronics Engineers (IEEE): Piscataway, NJ, USA, 2019; pp. 266–269. [Google Scholar]
Serpelloni, E.; Faccenna, C.; Spada, G.; Dong, D.; Williams, S. Vertical GPS ground motion rates in the Euro-Mediterranean region: New evidence of velocity gradients at different spatial scales along the Nubia-Eurasia plate boundary. J. Geophys. Res. Solid Earth 2013, 118, 6003–6024. [Google Scholar] [CrossRef] [Green Version]
Wegmuller, U.; Strozzi, T.; Bitelli, G. Validation of ERS differential SAR interferometry for land subsidence mapping: The Bologna case study. In Proceedings of the IEEE 1999 International Geoscience and Remote Sensing Symposium, IGARSS’99 (Cat. No.99CH36293), Hamburg, Germany, 28 June–2 July 1999; pp. 1131–1133. [Google Scholar] [CrossRef]
Solari, L.; Del Soldato, M.; Bianchini, S.; Ciampalini, A.; Ezquerro, P.; Montalti, R.; Raspini, F.; Moretti, S. From ERS 1/2 to Sentinel-1: Subsidence Monitoring in Italy in the Last Two Decades. Front. Earth Sci. 2018, 6, 149. [Google Scholar] [CrossRef]
Del Soldato, M.; Farolfi, G.; Rosi, A.; Raspini, F.; Casagli, N. Subsidence Evolution of the Firenze–Prato–Pistoia Plain (Central Italy) Combining PSI and GNSS Data. Remote Sens. 2018, 10, 1146. [Google Scholar] [CrossRef] [Green Version]
Gambolati, G. CENAS: Coastal Evolution of the Upper Adriatic Sea due to Sea Level Rise and Natural and Anthropogenic Land Subsidence; Kluwer Academic Publishers: Dordrecht, The Netherlands, 1998. [Google Scholar]
Lavecchia, G.; Castaldo, R.; De Nardis, R.; De Novellis, V.; Ferrarini, F.; Pepe, S.; Brozzetti, F.; Solaro, G.; Cirillo, D.; Bonano, M.; et al. Ground deformation and source geometry of the 24 August 2016 Amatrice earthquake (Central Italy) investigated through analytical and numerical modeling of DInSAR measurements and structural-geological data. Geophys. Res. Lett. 2016, 43, 12389–12398. [Google Scholar] [CrossRef]
Cheloni, D.; De Novellis, V.; Albano, M.; Antonioli, A.; Anzidei, M.; Atzori, S.; Avallone, A.; Bignami, C.; Bonano, M.; Calcaterra, S.; et al. Geodetic model of the 2016 Central Italy earthquake sequence inferred from InSAR and GPS data. Geophys. Res. Lett. 2017, 44, 6778–6787. [Google Scholar] [CrossRef]
Rutter, C.M.; Isaaks, E.H.; Srivastava, R.M. An Introduction to Applied Geostatistics. J. Am. Stat. Assoc. 1991, 86, 548. [Google Scholar] [CrossRef]
Fornaro, G.; Franceschetti, G.; Lanari, R.; Sansosti, E.; Tesauro, M. Global and local phase-unwrapping techniques: A comparison. J. Opt. Soc. Am. A 1997, 14, 2702–2708. [Google Scholar] [CrossRef]
Fornaro, G.; Franceschetti, G.; Lanari, R.; Rossi, D.; Tesauro, M. Interferometric SAR phase unwrapping using the finite element method. IEE Proc.—Radar Sonar Navig. 1997, 144, 266. [Google Scholar] [CrossRef]
Trasatti, E.; Casu, F.; Giunchi, C.; Pepe, S.; Solaro, G.; Tagliaventi, S.; Berardino, P.; Manzo, M.; Pepe, A.; Ricciardi, G.P.; et al. The 2004–2006 uplift episode at Campi Flegrei caldera (Italy): Constraints from SBAS-DInSAR ENVISAT data and Bayesian source inference. Geophys. Res. Lett. 2008, 35. [Google Scholar] [CrossRef] [Green Version]
Manconi, A.; Walter, T.R.; Manzo, M.; Zeni, G.; Tizzani, P.; Sansosti, E.; Lanari, R. On the effects of 3-D mechanical heterogeneities at Campi Flegrei caldera, southern Italy. J. Geophys. Res. Space Phys. 2010, 115. [Google Scholar] [CrossRef] [Green Version]
D’Auria, L.; Giudicepietro, F.; Martini, M.; Lanari, R. The 4D imaging of the source of ground deformation at Campi Flegrei caldera (southern Italy). J. Geophys. Res. Space Phys. 2012, 117. [Google Scholar] [CrossRef]
Samsonov, S.; Tiampo, K.F.; Camacho, A.G.; Fernández, J.; Gonzalez, P.J. Spatiotemporal analysis and interpretation of 1993-2013 ground deformation at Campi Flegrei, Italy, observed by advanced DInSAR. Geophys. Res. Lett. 2014, 41, 6101–6108. [Google Scholar] [CrossRef] [Green Version]
D’Auria, L.; Pepe, S.; Castaldo, R.; Giudicepietro, F.; Macedonio, G.; Ricciolino, P.; Tizzani, P.; Casu, F.; Lanari, R.; Manzo, M.; et al. Magma injection beneath the urban area of Naples: A new mechanism for the 2012–2013 volcanic unrest at Campi Flegrei caldera. Sci. Rep. 2015, 5, 13100. [Google Scholar] [CrossRef] [Green Version]
Ganas, A.; Elias, P.; Bozionelos, G.; Papathanassiou, G.; Avallone, A.; Papastergios, A.; Valkaniotis, S.; Parcharidis, I.; Briole, P. Coseismic deformation, field observations and seismic fault of the 17 November 2015 M = 6.5, Lefkada Island, Greece earthquake. Tectonophys. 2016, 687, 210–222. [Google Scholar] [CrossRef]
Bie, L.; Gonzalez, P.J.; Rietbrock, A. Slip distribution of the 2015 Lefkada earthquake and its implications for fault segmentation. Geophys. J. Int. 2017, 210, 420–427. [Google Scholar] [CrossRef]
Evstatiev, D.; Evlogiev, Y. Landslides along the northern Black Sea coast between Varna city and Kavarna town (Bulgaria). Geo-Eco-Marina 2013, 19, 39–58. [Google Scholar]
Tang, W.; Motagh, M.; Zhan, W. Monitoring active open-pit mine stability in the Rhenish coalfields of Germany using a coherence-based SBAS method. Int. J. Appl. Earth Obs. Geoinf. 2020, 93, 102217. [Google Scholar] [CrossRef]
Perski, Z. ERS InSAR data for Geological Interpretation of Mining Subsidence in Upper Silesian Coal Basin in Poland. In Proceedings of the FRINGE’99 2nd International Workshop on SAR Interferometry, Liege, Belgium, 10–12 November 1999. [Google Scholar]
Liu, X.; Hu, J.; Sun, Q.; Li, Z.W.; Zhu, J. Deriving 3-D Time-Series Ground Deformations Induced by Underground Fluid Flows with InSAR: Case Study of Sebei Gas Fields, China. Remote Sens. 2017, 9, 1129. [Google Scholar] [CrossRef] [Green Version]
Hoffmann, J.; Zebker, H.A.; Galloway, D.L.; Amelung, F. Seasonal subsidence and rebound in Las Vegas Valley, Nevada, observed by Synthetic Aperture Radar Interferometry. Water Resour. Res. 2001, 37, 1551–1566. [Google Scholar] [CrossRef]
Schreier, G. Opportunities by the Copernicus Program for Archaeological Research and World Heritage Site Conservation. Remote Sens. Archaeol. Cult. Landsc. 2020. [Google Scholar] [CrossRef]
Bonano, M.; Manunta, M.; Zinno, I.; Buonanno, S.; Ojha, C.; Lanari, R. A large scale exploitation of high resolution satellite SAR data to analyze surface deformation in urban areas through the parallel full resolution SBAS-DInSAR approach. In Proceedings of the 21st EGU General Assembly, Vienna, Austria, 7–12 April 2019. [Google Scholar]
Tarek, M.; Brissette, F.P.; Arsenault, R. Evaluation of the ERA5 reanalysis as a potential reference dataset for hydrological modeling over North-America. Hydrol. Earth Syst. Sci. 2020, 24, 2527–2544. [Google Scholar]
Zhu, X.X.; Tuia, D.; Mou, L.; Xia, G.-S.; Zhang, L.; Xu, F.; Fraundorfer, F. Deep Learning in Remote Sensing: A Comprehensive Review and List of Resources. IEEE Geosci. Remote Sens. Mag. 2017, 5, 8–36. [Google Scholar] [CrossRef] [Green Version]
Ma, L.; Liu, Y.; Zhang, X.; Ye, Y.; Yin, G.; Johnson, B.A. Deep learning in remote sensing applications: A meta-analysis and review. ISPRS J. Photogramm. Remote Sens. 2019, 152, 166–177. [Google Scholar] [CrossRef]
The European Ground Motion Service (EGMS). Available online: https://land.copernicus.eu/user-corner/technical-library/european-ground-motion-service (accessed on 10 September 2020).

Figure 1. Block diagram of the S-1 IWS P-SBAS processing chain workflow modified from [20].

Figure 2. Workflow of the automatic CC-based pipeline implementing the S-1 IWS P-SBAS processing. Green, blue and red arrows represent information exchange, data transfer and results visualization, respectively.

Figure 3. Block diagram of the proposed procedure for the DInSAR deformation signals refinement.

Figure 4. Examples of the GNSS stations position time series screening operation. Plots (1)–(5) are relevant to the five GNSS stations named ASCC, CESI, ALDB, BARY and GINA, respectively. On the left side the original time series (East-West component) are shown, while on the right side there are the corresponding ones following the data screening procedure. The dashed blue vertical lines identify the dates for which the GNSS auxiliary files report the presence of gaps and/or steps in the time series.

Figure 5. Planar components of the GNSS stations identified through our data screening and selection procedure, superimposed on the 1-arcsec STRM DEM of the investigated area. (A,B) The East/West velocity components before (A) and after (B) the implemented procedure. (C,D) The same as (A,B) but for the North/South velocity components.

Figure 6. Pictorial example of the implemented DInSAR deformation results refinement procedure applied to an S-1 IWS frame relevant to a portion of the Central-Northern Italy. (A) Mean deformation velocity map following the straightforward application of the P-SBAS approach. (B) Mean deformation velocity map obtained by interpolating the data relevant to the GNSS stations selected following the lines of Section 3.1. (C) Mean deformation velocity map of the high spatial frequency component. (D) Mean deformation showing the velocity map of the overall signals. The white star in panels A and B identifies the reference pixel, while the white arrows in the lower right corners show the SAR sensor azimuth/range directions. At the bottom a block diagram of the implemented procedure is shown.

Figure 7. SRTM DEM map of the investigated portion of the Central-Northern Italy, shown in Figure 6, highlighting the locations of the GNSS stations. Blue and yellow markers represent stations exploited for and excluded from the DInSAR deformation signals refinement step, respectively.

Figure 8. (A–R) Comparison between the P-SBAS (black triangles) and the LOS-projected GNSS (red stars) deformation time series relevant to the stations labelled in Figure 7 with capital letters from A to R. The standard deviation values (in cm) of the difference between the two time series are also reported.

Figure 9. SRTM DEM of the investigated area with the locations of the analyzed S-1 IWS frames. Note that the red and blue colors identify the frames whose images have been collected during the March 2015–September 2018 and March 2015–March 2020 time intervals, respectively.

Figure 10. Overall GNSS stations analysed in the presented study. Blue and yellow markers represent the exploited and the discarded stations, respectively, following the screening and selection procedure discussed in Section 3.

Figure 11. Processing architecture implemented within the ONDA DIAS CC Platform.

Figure 12. Mosaicking of the LOS mean deformation velocity maps (cm/year), geocoded and superimposed on an optical image of Europe (Mapbox Satellite Streets source). The white rectangles identify the zoom-in areas that are analyzed in more details in the following. Note that the red color corresponds to a sensor-target range increase, while the blue one to a sensor-target range decrease. The inset in the upper right corner reports the SAR sensor azimuth/range and the North/East directions.

Figure 13. Zoomed-in view of the mean deformation velocity map relevant to the Campi Flegrei caldera (Italy). (a) LOS mean deformation velocity map superimposed on the SRTM DEM of the zone. (b) Displacement time series relevant to the pixel P1, marked by the white star in the map, located in the maximum deforming area. (c) Plot of the mean deformation velocity along the section crossing the maximum uplift zone and identified by the white dashed line in panel (a).

Figure 14. Zoomed-in view of the mean deformation velocity map relevant to the Lefkada Island (Greece). (a) LOS mean deformation velocity map superimposed on the SRTM DEM of the zone. (b,c) Displacement time series relevant to the pixels P1 and P2, marked by the white stars in (a) and located in the maximum deforming co-seismic areas. The vertical dashed red lines identify the 17 November 2015 earthquake.

Figure 15. Zoomed-in view of the mean deformation velocity map relevant to the Northern Black Sea coast (Bulgaria). (a) LOS mean deformation velocity map superimposed on the SRTM DEM of the zone. (b,c) Displacement time series relevant to the pixels, marked by the white stars in (a) and located in the areas of Balchik town (P1) and Kranevo village (P2), both affected by landslides movements.

Figure 16. Zoomed-in views of the mean deformation velocity map relevant to two wide areas located in Germany (a) and Poland (b), both affected by extensive mining activities. (a,b) LOS mean deformation velocity map superimposed on the SRTM DEM of the zones; the white rectangles identify two areas that are shown in the panels (c,d). (c) Zoom of the LOS mean deformation velocity map in (a) relevant to the Nocthen coalmine (Germany). (d) Zoom of the LOS mean deformation velocity map in (b) relevant to the Katowice coalmine (Poland). (e) Displacement time series relevant to the pixel P1 marked by the white star in (c). (f) Displacement time series relevant to the pixel P2 marked by the white star in (d).

Figure 17. Correlation map between DInSAR time series and an annual sinusoid relevant to an area located in Northern Italy, close to the city of Milano. (a) Peak-to-peak oscillation amplitude map superimposed on the SRTM DEM of the zone. Note that pixels with peak-to-peak amplitude values equal to 0 are those for which the correlation coefficient is smaller than 0.8. (b–d) Displacement time series relevant to the pixels P1, P2 and P3, marked by the black stars in the map, located in correspondence to the Settala, Sergnano and Ripalta Guerina gas storage sites, respectively.

Figure 18. Correlation map between the DInSAR time series and an annual sinusoid relevant to the Firenze-Prato-Pistoia basin (North Italy). (a) Peak-to-peak oscillation amplitude map superimposed on the SRTM DEM of the zone. Note that pixels with peak-to-peak amplitude values equal to 0 are those for which the correlation coefficient is smaller than 0.8. (b) Displacement time series relevant to the pixel P1, marked by the black star in the map, located in correspondence of the maximum deforming zone.

Table 1. ONDA DIAS Computing Nodes Configuration.

ONDA DIAS Computing Nodes Configuration
Architecture	64 bit
vCPU/vCores *	64
RAM	512 GB
Internal Disks	18.5 TB (5 SSD in RAID 0 configuration)
Network	2 Gb/s

* vCPU and vCores are acronyms of virtual CPU and virtual Cores, respectively.

© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

Share and Cite

MDPI and ACS Style

Lanari, R.; Bonano, M.; Casu, F.; Luca, C.D.; Manunta, M.; Manzo, M.; Onorato, G.; Zinno, I. Automatic Generation of Sentinel-1 Continental Scale DInSAR Deformation Time Series through an Extended P-SBAS Processing Pipeline in a Cloud Computing Environment. Remote Sens. 2020, 12, 2961. https://doi.org/10.3390/rs12182961

AMA Style

Lanari R, Bonano M, Casu F, Luca CD, Manunta M, Manzo M, Onorato G, Zinno I. Automatic Generation of Sentinel-1 Continental Scale DInSAR Deformation Time Series through an Extended P-SBAS Processing Pipeline in a Cloud Computing Environment. Remote Sensing. 2020; 12(18):2961. https://doi.org/10.3390/rs12182961

Chicago/Turabian Style

Lanari, Riccardo, Manuela Bonano, Francesco Casu, Claudio De Luca, Michele Manunta, Mariarosaria Manzo, Giovanni Onorato, and Ivana Zinno. 2020. "Automatic Generation of Sentinel-1 Continental Scale DInSAR Deformation Time Series through an Extended P-SBAS Processing Pipeline in a Cloud Computing Environment" Remote Sensing 12, no. 18: 2961. https://doi.org/10.3390/rs12182961

APA Style

Lanari, R., Bonano, M., Casu, F., Luca, C. D., Manunta, M., Manzo, M., Onorato, G., & Zinno, I. (2020). Automatic Generation of Sentinel-1 Continental Scale DInSAR Deformation Time Series through an Extended P-SBAS Processing Pipeline in a Cloud Computing Environment. Remote Sensing, 12(18), 2961. https://doi.org/10.3390/rs12182961

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Menu