Decomposing DInSAR Time-Series into 3-D in Combination with GPS in the Case of Low Strain Rates: An Application to the Hyblean Plateau, Sicily, Italy

Abstract

Differential Interferometric SAR (DInSAR) time-series techniques can be used to derive surface displacement rates with accuracies of 1 mm/year, by measuring the one-dimensional distance change between a satellite and the surface over time. However, the slanted direction of the measurements complicates interpretation of the signal, especially in regions that are subject to multiple deformation processes. The Simultaneous and Integrated Strain Tensor Estimation from Geodetic and Satellite Deformation Measurements (SISTEM) algorithm enables decomposition into a three-dimensional velocity field through joint inversion with GNSS measurements, but has never been applied to interseismic deformation where strain rates are low. Here, we apply SISTEM for the first time to detect tectonic deformation on the Hyblean Foreland Plateau in South-East Sicily. In order to increase the signal-to-noise ratio of the DInSAR data beforehand, we reduce atmospheric InSAR noise using a weather model and combine it with a multi-directional spatial filtering technique. The resultant three-dimensional velocity field allows identification of anthropogenic, as well as tectonic deformation, with sub-centimeter accuracies in areas of sufficient GPS coverage. Our enhanced method allows for a more detailed view of ongoing deformation processes as compared to the single use of either GNSS or DInSAR only and thus is suited to improve assessments of regional seismic hazard.

1. Introduction

The use of space-borne Synthetic Aperture Radar (SAR) imagery has had considerable success during the last two decades for the monitoring of the seismic cycle [1,2]. Due to the coherent measurement principle of a SAR system, the interference of the signal of two SAR images (i.e., Interferometric SAR-InSAR) can be calculated. Using a reference of the a priori distance between the satellite and the sensed object in form of a Digital Elevation Model (DEM), the differential phase of 2 images, acquired at different dates, can be exploited (i.e., Differential InSAR-DInSAR). This results in measurements of change in distance between the antenna and the ground in the 1-dimensional line-of-sight (LOS) direction of the radar aperture, with centimeter accuracies for the time span of the 2 acquisitions [3].

More recent developments of time-series DInSAR methods exploit the temporal evolution of the phase signal in time by using a stack of multiple SAR images, whereby the respective accuracies increase significantly. Advanced techniques referred to as Persistent Scatterers (PS) [4,5,6] and Small BASeline (SBAS) [7,8,9] have already demonstrated their potential in assessing ground surface motion for active tectonic zones such as plate boundaries (e.g., [10,11]), volcanoes (e.g., [12,13,14]), landslides (e.g., [15]) as well as urban areas (e.g., [16,17,18,19,20,21]), where local deformation rates exceed the expected error and enough pixels per area contain a relatively stable phase in time (i.e., coherent areas). Instead, the application of time-series DInSAR techniques for the detection of local seismicity at secondary faults of wider plate boundary areas is more challenging, since respective ground displacements are usually less than 1 cm/year. Under ideal circumstances, accuracies of 1 mm/year and less are reported for LOS measurements using time-series DInSAR methods [19,21,22]. This, however, is often not the case.

Furthermore, detected ground displacements in LOS are directional ambiguous and can originate both from horizontal as well as vertical deformations, which, in turn, may be results of different processes taking place. This hampers the 3-dimensional characterization of tectonic processes and may mislead their differentiation from other drivers of surface deformation such as anthropogenic induced subsidence due to construction, mining or water pumping.

In order to enhance robustness and enable the use in operational scenarios of (seismic) hazard assessments, two challenges need to be overcome: (1) to account properly for atmospheric noise due to changes in the atmospheric state between two or more SAR acquisitions and (2) to resolve the ambiguity of the measurements with respect to the direction of the ground movement due to the 1-dimensionality of the LOS measurement.

In this paper, we present a combination of recent advances in order to address both challenges over the area of the Hyblean plateau in south-east Sicily, Italy. This region experienced the largest earthquake (Mw $~$7.4) that struck Italy during historical times. The mainshock on the 11th of January 1693 was accompanied by a large tsunami wave affecting the eastern Sicilian coast [23]. As the north-easternmost part of the Pelagian block the area is located in the Eurasian-Nubian convergence zone along the Maghrebian belt and is characterized by complex geodynamic processes. Additionally, karst processing on the carbonatic massif [24] and water pumping for agricultural and industrial water supply may deform the surface locally [25]. While local studies of time-series DInSAR already showed anthropogenic induced subsidence along the coast [25], there have not been any analyses by means of combined usage of Global Positioning System (GPS) and time-series DInSAR over the whole region so far.

1.1. Techniques of Atmospheric Correction

Atmospheric noise forms one of the largest challenges faced by the InSAR community as it typically masks smaller magnitude signals [2,3,26,27,28]. The atmospheric noise is split in an ionospheric and tropospheric component. Free electrons in the ionosphere lead to partial dispersion of the electro-magnetic wave, which generates an apparent noise signal within an interferogram. This influence is however frequency-dependent and becomes more prominent in the longer wavelength spectra of L- and P-Band. For our study, which deals with standard image mode C-band imagery from the ENVISAT Advanced SAR (ASAR) instrument, the effect is expected to be small and we neglect it. It should be noted, however, that for recently launched C-band sensors covering a wider area (e.g., Sentinel-1), those effects are more significant.

Tropospheric path delays, instead, are caused by differences in temperature, atmospheric pressure, and foremost the water vapor (i.e., relative humidity). Those differences introduce a noise signal on various spatial scales [29] and can exhibit up to a 14 cm error in a single interferogram [3]. As tropospheric delays are integrated along the radar traveled path, the introduced noise can exhibit a strong topography-correlated component as well. Different methods exist that either estimate the tropospheric correction based on the interferometric phase or by using auxiliary data.

Traditionally, phase-based corrections are used by applying a spatiotemporal filter on DInSAR stacks. This approach aims to reduce the spatially-correlated atmospheric noise, assuming that the atmosphere behaves randomly in time while being correlated in space [6,9,14]. Other attempts estimate linear and power-law relationships between the interferometric phase and the topography, with the aim to reduce the topography-correlated as well as the long-wavelength tropospheric noise [29,30]. For all of those methods the challenge lies in avoiding leakage of the deformation signal in the atmospheric correction. Alternative approaches based on auxiliary data include corrections estimated from GPS data [31,32], spectrometer measurements [33], meteorological model data [27,34,35] as well as combined approaches [36,37,38]. While each of these techniques is capable of reducing the tropospheric InSAR noise, the success rate is often limited by the lower spatial resolution or availability of the auxiliary data [29].

Here we utilize the Toolbox for Reducing Atmospheric InSAR Noise (TRAIN, www.davidbekaert.com), which includes all the current state of the art methods [39]. For our study we chose a correction technique based on the 75 km resolution ERA-INTERIM weather model data from the European Center for Meteorological Weather Forecast (ECMWF) [40]. Due to the coarse resolution, weather models often fail to correctly represent the turbulent atmosphere, but they do reduce the topography-correlated and long-wavelength signal on average [39,41,42]. Therefore, diurnal, semi-diurnal as well as seasonal signals can be reduced within single interferograms. Furthermore, the approach holds the advantage of being non-parametric and thus can be applied fully automatic. Since the dataset provides global atmospheric reanalysis from 1979 onwards, it is universally applicable to all past, present and future InSAR data, hence providing improved estimates of the LOS measurements.

1.2. Techniques for Decomposing the DInSAR Signal into a 3-Dimensional Velocity Field

Regional seismic hazard assessments greatly benefit from a decomposed (i.e., 3-dimensional) velocity field as provided by geodetic GPS measurements due to the directional unambiguity in the assessment of ground motion. In addition, it becomes easier to distinguish tectonic processes from other causes of surface deformations. While DInSAR provides a much higher spatial resolution as a GPS network, it is limited to the 1-dimensionality of the LOS measurement, being a superposition from horizontal and vertical displacements.

Decomposing multiple tracks of DInSAR observations is not trivial and has been subject to recent research activities [43]. Once deformation rates are decomposed in the horizontal and vertical directions, parameters like slip rate, dilatation or strain can be estimated more precisely. Those parameters, in turn, are of importance for defining the boundary conditions of the respective faults within seismic hazard models [44].

In general, a first aspect takes into consideration the different viewing geometries of the descending and ascending orbits, which allows for untangling the horizontal from the vertical motion. However, due to the orbital configuration of space-borne SAR systems, the signal remains almost insensitive to north-south directional movements. By neglecting north-south displacements, it becomes mathematically feasible to disentangle the east-west and up-down directions from DInSAR data alone [45].

Pixel offset tracking [46] and Multiple-Aperture Interferometry (MAI) [47,48] are alternative SAR based techniques that allow for displacement mapping along the azimuth direction of the satellite flight path and exhibit more sensitivity to north-south displacements (e.g., [49,50,51]). Accuracies are however on the order of a few centimeters, which in most cases exceed the actual deformation rate of interseismic motion.

The integration of evenly distributed GPS stations together with DInSAR datasets is considered the only appropriate method to obtain a spatially dense, 3-dimensional velocity field with sub-centimeter accuracies [43]. While GPS measurements usually provide precise measures with millimeter accuracies in all 3 directions of the velocity field, the distance between the GPS stations is often too large to capture local variations in the tectonic signal. The complementary use of both methods allows combining the 3-dimensional information of the GPS data with the dense spatial sampling of the DInSAR and thus to overcome the limitations of both techniques.

While various methods have been proposed to integrate the two different observations, there is no consensus yet on an ideal strategy [52]. The main issue relates to the question of how surrounding GPS measurements are used to disentangle the LOS displacements for pixels that do not overlap with the GPS stations. First attempts were based on an a priori interpolation of the GPS velocity field. Its consequent integration with the DInSAR datasets from ascending and descending orbits is based on either a Bayesian statistical approach in combination with a Markov random field [53] or an analytical optimization [11]. Both approaches have in common that the optimization procedures used are not based on the physics of the deformation. Additionally, both rely on the result of the a priori interpolation usually done by means of the kriging technique for each component. This technique, in turn, requires an appropriate choice of a theoretical semivariogram model for each component which depicts one of the main critics in geostatistics [54]. A slightly different method has been proposed by Catalão et al., 2011, which uses the horizontal displacements derived from interpolated GPS measurements to subsequently calculate the vertical displacements from the time-series DInSAR in ascending and descending orbits separately. The latter are then adjusted to the vertical GPS velocities by minimizing the bias between both observations and an improved estimate of the vertical velocity field is obtained [55].

Wang & Wright [56] adapted a method that was previously developed for GPS measurements solely [57]. The original method is based on a triangulation network, whereas a constant linear motion gradient in space is assumed. The data integration is based on the interpolation of the GPS and DInSAR data by using a biconjugate gradients method. Additionally, long-wavelength orbital and topographical correlated atmospheric noise is attributed and a Laplacian smoothing operator is incorporated to remove short-wavelength features of non-tectonic motion. While the method shows great potential for large-scale applications with regard to tectonic motion only, the velocity field model is calculated for the vertices of the triangles. Thus, the resolution of the original InSAR data is reduced to the distance of the vertices and small-scale phenomena in fragmented areas like the Hyblean foreland might be smoothed out.

For our study we have tested the Simultaneous and Integrated Strain Tensor Estimation from Geodetic and Satellite Deformation Measurements (SISTEM) approach [58], which is capable of simultaneously integrating various geodetic measurements to disentangle the 1-dimensional LOS displacements estimates and derive a spatially dense, 3-dimensional velocity field. With respect to the other techniques it features the advantage that it operates on the full resolution of the DInSAR data while no preliminary interpolation of the GPS data is necessary. Furthermore, it follows the physically based concept of the elastic theory, which is commonly applied for interseismic strain accumulation [59] as well as aquifer compression where intergranular pressure is increased due to water pumping [60]. As demonstrated by Guglielmino et al. [61] the integration of additional geodetic measurements such as leveling and DInSAR time-series from different sensors is straightforward. Hence, one can take full advantage of available data sources over a certain area of interest. SISTEM has already been applied successfully to volcanic [58] as well as for coseismic events [61,62], but has not before been applied to the derivation of small deformation rates associated with interseismic strain accumulation. This study aims to demonstrate its usefulness also for such types of deformation by concurrently providing information on other types of displacements such as anthropogenic induced subsidence.

It has been shown that under the assumption of infinitesimal and homogeneous strain, geodetic benchmark measurements can be used to derive the full strain field adopting a linear approach [63,64]. Under this assumption the 3-dimensional displacement components $U_i(i=1\dots 3)$ of an arbitrary point P with the position $x_0=(x_{10},x_{20},x_{30})$ can be estimated by using N surrounding displacements measurements $u_{\left(n\right)}=(u_{1\left(n\right)},u_{2\left(n\right)},u_{3\left(n\right)})$ of the experimental points ($EP$) with the positions $x_n=(x_{1\left(n\right)},x_{2\left(n\right)},x_{3\left(n\right)})$ in the form of N equations:

$$u_{i\left(\right)n}=H_{ij}{\Delta }^{ij}{x}_{j\left(n\right)}+U_i(i,j=1\dots 3)$$

(1)

where $\delta x_{j\left(n\right)}=x_{j\left(n\right)}-x_{j\left(0\right)}$ describes the distance between the nth $EP$ and P, while $H_{ij}=\delta u_i/\delta x_i$ represents the elements of the displacement gradient tensor which can be broken down into the symmetric (i.e., the strain tensor E) and asymmetric part (i.e., the rigid body rotation tensor Ω) in the form of $H=E+\Omega$.

For every arbitrary point, $3N$ equations of type (1) can be transformed into the form of a linear matrix equation $A\times l=u+e$ and subsequently solved by the weighted least-square (WLS) approach. A represents the $3N\times 12$ design matrix, l refers to the column vector of unknown parameters consisting of the 3D displacements and the elements of the strain tensor as well as the rigid body rotation tensor and u is the column observation vector with the displacement measurements of the N geodetic benchmarks. The term e is the residual vector that contains the standard errors of the solution due to the stochastic nature of the estimation problem.

It should be noted that the inverse of the covariance matrix C for the WLS solution is modified by introducing a weighting factor as described by Shen et al. [64] in the form of:

$$C^{\prime }=Cexp\left(-\frac{d_{\left(n\right)}}{d_0}\right)$$

(2)

where $d_{\left(n\right)}$ represents the distance between the nth $EP$ and the arbitrary Point P. In this way, the nearby $EP$s have a higher influence on the estimation of P.

Instead of the weighted interpolation between the geodetic benchmarks, SISTEM integrates the DInSAR measurements into the framework of the linear matrix equations by relating a number K of DInSAR datasets to the components of the displacements vector of an arbitrary point P as follows:

$$D_{LOS}^{SP}=[U_1,U_2,U_3]{\left[\begin{array}{c}{S}_x^{SP},S_y^{SP},S_z^{SP}\end{array}\right]}^T(S=1\dots K)$$

(3)

where $D_{LOS}^{SP}$ is the Sth known LOS displacement of the point P and $[S_x^{SP},S_y^{SP},S_z^{SP}]$ is the unit vector pointing from point P towards the Sth satellite acquisition. Consequently, the final design matrix A is constructed by $3N+K$ rows and has the following form:

$$A=\left[\begin{array}{cccccccccccc}1& 0& 0& \Delta x_{1\left(1\right)}& \Delta x_{2\left(1\right)}& \Delta x_{3\left(1\right)}& 0& 0& 0& 0& \Delta x_{3\left(1\right)}& -\Delta x_{2\left(1\right)}\\ 0& 1& 0& 0& \Delta x_{1\left(1\right)}& 0& \Delta x_{2\left(1\right)}& \Delta x_{3\left(1\right)}& 0& -\Delta x_{3\left(1\right)}& 0& \Delta x_{1\left(1\right)}\\ 0& 0& 1& 0& 0& \Delta x_{1\left(1\right)}& 0& \Delta x_{2\left(1\right)}& \Delta x_{3\left(1\right)}& \Delta x_{2\left(1\right)}& -\Delta x_{1\left(1\right)}& 0\\ .& .& .& .& .& .& .& .& .& .& .& .\\ .& .& .& .& .& .& .& .& .& .& .& .\\ .& .& .& .& .& .& .& .& .& .& .& .\\ 1& 0& 0& \Delta x_{1\left(N\right)}& \Delta x_{2\left(N\right)}& \Delta x_{3\left(N\right)}& 0& 0& 0& 0& \Delta x_{3\left(N\right)}& -\Delta x_{2\left(N\right)}\\ 0& 1& 0& 0& \Delta x_{1\left(N\right)}& 0& \Delta x_{2\left(N\right)}& \Delta x_{3\left(N\right)}& 0& -\Delta x_{3\left(N\right)}& 0& \Delta x_{1\left(N\right)}\\ 0& 0& 1& 0& 0& \Delta x_{1\left(N\right)}& 0& \Delta x_{2\left(N\right)}& \Delta x_{3\left(N\right)}& \Delta x_{2\left(N\right)}& -\Delta x_{1\left(N\right)}& 0\\ S_x^{\left(1\right)P}& S_y^{\left(2\right)P}& S_z^{\left(3\right)P}& 0& 0& 0& 0& 0& 0& 0& 0& 0\\ S_x^{\left(K\right)P}& S_y^{\left(K\right)P}& S_z^{\left(K\right)P}& 0& 0& 0& 0& 0& 0& 0& 0& 0\end{array}\right]$$

(4)

It can be noted that the first $3N$ rows refer to the N blocks of three equations containing the positional information of the 3D geodetic measurements by GPS. The last K rows, instead, correspond to the DInSAR data. Consequently, the observation vector u and the column vector of unknown parametersl are defined as follows:

$$u={[u_{1\left(1\right)},u_{2\left(1\right)},u_{3\left(1\right)}\cdots u_{1\left(n\right)},u_{2\left(n\right)},u_{3\left(n\right)},D_{LOS}^{1P},D_{LOS}^{KP}]}^T$$

(5)

$$l={[U_1,U_2,U_3,{ϵ}_{11},{ϵ}_{12},{ϵ}_{13},{ϵ}_{21},{ϵ}_{22},{ϵ}_{23},{ϵ}_{31},{ϵ}_{32},{ϵ}_{33},{\omega }_1,{\omega }_2,{\omega }_3]}^T$$

(6)

where ${ϵ}_{ij}$ represents the components of the strain tensor E defined as:

$$E=\left[\begin{array}{ccc}{ϵ}_{11}& {ϵ}_{21}& {ϵ}_{31}\\ {ϵ}_{12}& {ϵ}_{22}& {ϵ}_{32}\\ {ϵ}_{13}& {ϵ}_{23}& {ϵ}_{33}\end{array}\right]$$

(7)

and ${\omega }_i$ represents the components of the rigid body rotation tensor Ω defined as:

$$\Omega =\left[\begin{array}{ccc}0& -{\omega }_3& {\omega }_2\\ {\omega }_3& 0& -{\omega }_1\\ -{\omega }_2& {\omega }_1& 0\end{array}\right]$$

(8)

By combining SISTEM with the tropospheric corrected time-series DInSAR data, we demonstrate that sub-centimeter measurements of the 3-dimensional velocity field are feasible on a routine basis if sufficient spatial coverage of the geodetic input data is given. Such information can then be used to identify previously unknown surface deformation as well as characterizing and constraining seismic source parameters for respective hazard assessments in the presence of tectonic activity.

2. Method and Data Collection

2.1. Study Area

The Hyblean plateau is located in the south-east of the Italian island of Sicily, in the center of the Mediterranean Sea (Figure 1). The region features a complex tectonic setting of the Eurasia-Nubian convergence zone along the Maghrebian Belt. A manifold geodynamic evolution and a set of different crustal blocks trapped between the Eurasian and African rigid plates, led to a fragmented and irregularly shaped convergence zone of the African and European plate [65,66]. According to Ben Avraham & Grasso 1991 [67], the plateau is considered as a block of the African foreland that resists the underthrusting beneath the Maghrebian belt due to its high buoyancy. In turn, this formed the arc-shaped deflection of the Maghrebian belt along the northern edge of the Hyblean foreland.

On the eastern side, this block is bounded by the Hyblean-Maltese escarpment, a normal NNW-SSE trending and ENE dipping fault system that decouples the thick inland continental crust from the Ionian Mesozoic oceanic crust underthrusting the Calabrian arc [23]. The Ionian margin is thereby characterized by a typical horst and graben sequence, oblique with respect to the Hybleo-Maltese escarpment.

Downbending of the foreland takes place along the Gela-Catania foredeep. NE-SW trending normal faults characterize the area and form wide structures along the Scordia-Lentini graben. On the western side, the Scicli-Ragusa fault extends from the Scordia-Lentini basin in the north towards the Marina di Ragusa graben in the south and divides the foreland into two crustal blocks. The NE-SW normal faults controlling the Scordia-Lentini and the Marina di Ragusa extensional basins as well as the fault segments of the Scicli line have been partially reactivated by both reverse and left-lateral motions, respectively [70,71,72].

Historically, major destructive earthquakes hit the region in 1196, 1542 and 1693, though their location is still debated [23,73,74]. The last major earthquake (M${}_L$ = 5.4) occurred on 13 December 1990 located at the offshore part of the Lentini Graben [75,76]. The understanding of the neotectonics is of major interest for assessments regarding the regional seismic hazard. According to the national hazard maps it belongs to one of Italy’s most endangered regions with a 10% probability of exceeding a peak ground acceleration of more than 0.25 g within the next 50 years [77]. As shown in Figure 1b, there have been more than 450 small to medium strong earthquakes (M${}_L$ > 1) throughout the last 2 decades [66,78]. Despite numerous research activities it is still debated how the internal deformation of Sicily is accommodated by the faults and the number of faults that may take up this deformation [71].

2.2. Data Collection & Processing

2.2.1. DInSAR Processing

The multi-temporal DInSAR analysis has been performed using 49 ascending and 58 descending Envisat ASAR images that cover the time period from 2003 until 2010 (Figure 2). The SAR imagery has been focused using the ROI_PAC processor version 3.01 [79]. Subsequently, the interferograms were generated with the DORIS software [80]. We used the 30 m SRTM data for the topographic phase removal [81]. For the time-series processing the StaMPS Small Baseline (SB) approach was applied [7,82]. Due to the specific selection process of stable targets, this approach favors the detection of distributed scatterers which occur usually in rural regions. Furthermore, we could validate the phase unwrapping procedure by comparing the residuals of the phase between the SB interferograms and the modeled inversion for the single-master interferogram with the closed-loop test [7,8]. Noisy interferograms with spatially correlated phase unwrapping errors have been subsequently sorted out and the mean LOS velocities have been recalculated. In total we obtained a network of 119 interferograms for the descending and 148 interferograms for the ascending track. We then decreased further noise by resampling the pixel size to 300 m. The contribution related to associated DEM errors has been removed as part of the StaMPS time-series inversion. We corrected for the Envisat oscillator drift as described by Marinkovic & Larsen 2013 [83]. We further accounted for the long-wavelength orbit errors in our decomposition by using the GPS as constrain.

2.2.2. Tropospheric Correction

We corrected for tropospheric delays in the single interferograms by using ERA-INTERIM (75 km resolution) weather model data from the European Centre for Meteorological Weather Forecast (ECWMF) [40], which provides daily estimates of pressure, temperature and relative humidity at 0:00 h, 6:00 h, 12:00 h and 18:00 h UTC. These model outputs are used to estimate the tropospheric delay in the radar LOS for each interferogram. We used the routines available in the Toolbox for Reducing Atmospheric InSAR Noise (TRAIN) [84] described in Bekaert et al. 2015 [39]. The procedure performs a horizontal and vertical spline interpolation for the pressure, temperature and relative humidity. The slant tropospheric delay is then calculated by integrating the refractivity from the surface upwards, projected into the radar LOS. A linear interpolation is done in the temporal domain to match the SAR acquisition time of 22:00 h UTC for the ascending pass and 10:00 h UTC for the descending pass.

2.2.3. Spatial Enhancement

The coverage of the DInSAR datasets depends on the phase stability of every pixel. Since SISTEM only considers valid DInSAR pixels and does not interpolate per se over areas without information, a combined gap-filling and subsequent filtering procedure has been applied to enhance the spatial coverage of the DInSAR data prior to the data integration. This approach has been originally used to improve interferometric elevation data from the SRTM mission for land surface process applications [85]. The underlying problem is somewhat similar, since small-sized no-data areas are filled and residual artifacts from phase noise are removed by simultaneously preserving actual features like edges and ditches.

In our case, we assumed that for nearby areas the expected deformation signals are spatially correlated. As a first step, the gap-filling, which is based on a spline interpolation, has been applied for areas with less or equal than 100 pixels. In order to remove local artifacts that feature sort of a “salt and pepper” noise, we subsequently applied the multi-directional Lee filter [86]. The filter operates in the direction of the minimal variance across 16 directions within a windows size of 9 by 9 pixels. The actual distance of the filter operation is weighted by a parameter that takes the local noise into account.

2.2.4. GPS Data

GPS based ground deformation measurements over the Hyblean plateau have been routinely carried out since 1991 as a consequence of the M${}_L$ = 5.4 Augusta earthquake. The data used in this work has been collected during different geodetic surveys between 1993 and 2006 [70,87] as well as from continuously running GPS stations lying on the Hyblean plateau (Figure 3). We used the final precise velocities of 24 selected stations matching the DInSAR coverage. The velocities have been calculated and reported in Angelica et al. 2013 [65] and Mastrolembo Ventura et al. [71].

2.2.5. SISTEM Integration

SISTEM simultaneously decomposes the input measurements into a 3-dimensional field vector based on the elasticity theory. The original SISTEM approach, as described by Guglielmino et al. [58], considers only one DInSAR dataset that is combined with the GPS data. The simultaneous integration of GPS and DInSAR measurements into a linear matrix equation is then solved by the WLS method. However, SISTEM can be easily extended by other DInSAR datasets as well as other geodetic measurements. Since the algorithm calculates every pixel singularly, SISTEM is capable of adapting the DInSAR-GPS integration according to the number of available DInSAR measurements for a particular pixel. In our case, conditioned by the expected small displacement rates, preference was given to the highest possible constrain of the WLS solution considering only overlapping pixels of both DInSAR datasets used.

Another important factor is the weighting of the GPS measurements based on the distance-decay function with respect to the nearest stations, expressed as the level of locality. This level of locality drops after a certain distance threshold. Beyond this radius the DInSAR data becomes the dominant information source. The threshold can be set for every component of the 3-dimensional velocity vector separately. For this study, the thresholds have been set to 4, 2 and 1 km for the north-south, east-west and up-down component, respectively. The highest threshold for the north-south component has been chosen because it is theoretically not possible to derive this component from the DInSAR data due to the acquisition geometry. Information outside this radius, which is predominantly constrained by the DInSAR data, becomes therefore unreliable. The threshold of 2 km for the east-west component instead was chosen in order to pick up a higher degree of spatial variation by utilizing the DInSAR data. For the vertical component the lowest threshold was used for 2 reasons. Firstly, DInSAR is most sensitive to vertical displacements due to its steep look-angle (∼23 degree). Secondly, the GPS data exhibits the highest uncertainties in the vertical direction. This is especially true for the GPS stations that have been installed in the 1990’s for which the standard deviations within the vertical direction ranges from 1.2 to 3 mm/year so that the DInSAR measurements have been considered as equally accurate, if not better.

For the final computation the structure of the covariance matrix has to be modified by adding the variances of the observations. Instead of estimating them on the DInSAR and GPS measurements itself, fixed values have been used for computational effectiveness. In order to account for the spatially correlated noise, a fixed variance of 7 mm/year has been used for the DInSAR data, while for the horizontal and vertical GPS measurements the variances were set to 3 mm/year and 1 cm/year, respectively. This, however, entails biases in the estimations of the standard error of the solution. While areas of higher standard errors indicate either non-elastic displacements or poor input data quality, they are not suited to determine the actual accuracy of the output data due to the arbitrarily chosen variances. In order to evaluate the overall performance of the data integration, the 3-dimensional velocity field was therefore used to recalculate the LOS velocities. The resultant residuals with respect to the input LOS measurements give a proper indication of the actual accuracy.

3. Results & Discussion

3.1. Tropospheric Correction

In order to validate the expected improvements of the mean LOS velocities after the tropospheric correction, two measures have been carried out. First, we compared the GPS velocities projected to the LOS of both tracks with our DInSAR results. We find a reduction of the Root Mean Square Error (RMSE) that amounts to about 15% for both tracks (

Table 1. Accuracy of DInSAR results before and after tropospherical correction. RMSE of GPS velocities projected into the LOS against the DInSAR results is given (middle column) as well as the mean RMSE for all interferograms used in the LOS velocity inversion (right column).

	Mean LOS vs. GPS (RMSE)	Mean RMSE of all IFGs
non-corrected Stack (Ascending)	3.4 mm/year	3.9 rad
Tropos. corrected stack (Ascending)	2.8 mm/year	3.5 rad
Percentual Improvement	17.6%	10%
non-corrected Stack (Descending)	3.3 mm/year	4.2 rad
Tropos. corrected stack (Descending)	2.8 mm/year	3.8 rad
Percentual Improvement	15%	9.5%

). The high RMSE of 2.8 mm/year is conditioned by stations with older GPS instrumentation and larger uncertainties foremost in the vertical direction. By excluding these stations, the RMSE drops to 1.4 mm/year for the descending, and 1.5 mm/year for the ascending track, whereas the absolute residuals range from 0.1 mm/year up to 1.5 mm/year. Additional confirmation comes from the assessment of the variability of the mean LOS measurements within each single interferograms by calculating the RMSE. A reduction in RMSE of the tropospherically-corrected interferogram with respect to the uncorrected one indicates a general smoothing of the data, which would be expected due to the partial elimination of the tropospheric noise. Indeed, a reduction of about 10% could be achieved for the interferograms of both tracks on average (Table 1), indicating a decrease of long-wavelength tropospheric noise. These improvements are on the order of what was observed for a similar study over Eastern Turkey [41].

The generally smoother deformation field of the tropospherically corrected mean LOS velocities is also visible in the LOS displacement maps (Figure 4b,f). The eastern part of the Plateau, in particular, shows a decrease in variability, which further indicates the improvement.

It is worth mentioning that we also tested high-resolution NWM data as proposed in Catalão et al., 2011, but similar to the results of Bekaert et al., 2015 over central Italy, the higher resolution did not improve our estimates.

3.2. Spatial Enhancement

In contrast to prior studies that applied SISTEM, we found that when observations are not present for both, ascending and descending data, artifacts appear within our inversion result. In order to avoid a significant drop in spatial coverage when limiting the observations to the mutual coverage of both tracks, we performed a preceding gap-filling procedure. Under the assumption that the detected displacements and the DInSAR signal are spatially correlated, the gap-filling procedure led to a significant increase in the number of overlapping pixels of the two datasets (Figure 4c,g). Since the study area is structurally fragmented, a spatial area threshold of 100 pixels (∼9 km${}^2$) has been set as a maximum size for the gap closed areas.

However, isolated pixels within such enclosed areas often suffer from increased phase noise so that in some of the gap filled areas a “salt and pepper” noise becomes apparent due to the higher weight in the interpolation process of those pixels. In order to smooth those areas, the multi-directional Lee filter was applied subsequently (Figure 4d,h). Since it is an anisotropic filtering approach, gradients along faults have not been affected as it would be the case for standard mean filters. By integrating the combined application of the gap-filling and anisotropic filtering in our workflow, it became possible to exploit the spatial correlation of the DInSAR data and enhance the coverage for the subsequent SISTEM data integration without introducing further noise. Indeed, the difference of the RMSE between the GPS velocities projected in the LOS and the tropospherically-corrected DInSAR LOS before and after this procedure amounts to 0.02.

3.3. SISTEM Decomposition

By integrating the GPS velocities (Figure 3) with the corrected and spatially enhanced mean LOS velocities (Figure 4d,h) into the SISTEM framework, the complementary features of both measurements are utilized to retrieve three-dimensional displacement estimates in the north-south, east-west and up-down direction (Figure 5a–c). The corresponding standard errors from the WLS solution are also provided as well (Figure 5a’–c’). In addition, the absolute residuals with respect to the original GPS data have been computed (Figure 5). However, this information should not be regarded as an independent validation source since the GPS data were used as an input into the decomposition. It rather indicates how the GPS measurements are affected by the integration with the DInSAR dataset. The standard error, instead, is biased by the fixed covariance values used for the WLS solution. The estimates are therefore unreliable in absolute terms, but can indicate relative errors that are related either to data input issues, or the circumstance that the surface is not deforming in an elastic manner. In order to realistically access the overall performance of SISTEM in terms of numerical accuracy, the retrieved velocity field has been used to recalculate the LOS velocities of the two SAR tracks. The residuals between the predicted LOS velocities from the SISTEM estimates and the original LOS (Figure 6) reveal a mean value of −0.45 mm/year with a standard deviation of 1.12 mm/year for the ascending, and a mean value of 0.76 mm/year with a standard deviation of 1.67 mm/year for the descending track. This result, in addition to the aforementioned deviations, ranging from 0.1 mm/year to 1.5 mm/year, between the GPS and the original LOS velocities into account, demonstrate the ability of this method to retrieve three-dimensional sub-centimeter surface deformation.

3.3.1. Vertical Component

As can be seen from the vertical displacement field (Figure 5c), the covered area is characterised by various zones of subsidence and uplift ranging from −10 to 5 mm/year. The comparison of the SISTEM output with the GPS data reveals the highest deviations of all three components. As mentioned earlier, older GPS stations suffer from poor vertical estimates as in the case of the ICRC station on the Monte Tauro Horst, northeast of the Augusta basin. Despite a significant uplift of about 6 mm/year indicated by the GPS velocity, the subsidence estimate of SISTEM in the neighbouring Augusta basin predominantly resembles the LOS displacements. Due to the high sensitivity of the DInSAR data with respect to the vertical component, and the lower weight given by the distance-decay function to the GPS station, SISTEM reinforces the respective estimation and reduces the magnitude of the standard error (Figure 5c’), as shown as well in the original approach [58]. Both tracks exhibit a LOS displacement of down to −10 mm/year for the Augusta basin, which translated into pure vertical movement, would result in about −11 mm/year for the acquisition geometry of ENVISAT Image mode. It is therefore underlined that despite the vicinity of a GPS station with contradictory vertical estimates, the SISTEM output is only marginally influenced by the GPS measurement, and thus improves the vertical displacement estimates using more accurate information provided by the mean LOS velocities of both DInSAR datasets. Indeed, this, and other areas affected by uplift and subsidence reveal low residuals of −1.5 to 1.5 mm/year when comparing the original LOS velocities to the predicted ones from the SISTEM decomposition (Figure 6). Nonetheless, GPS stations affected by low vertical accuracy have been kept for to their valuable, and more accurate horizontal components.

3.3.2. Horizontal Components

The retrieval of the north-south component (Figure 5a) as well as the east-west component (Figure 5b) is more critical since the DInSAR is less sensitive to this component. A good coverage of GPS stations is crucial for the retrieval of reliable estimates, especially for the north-south component, where the DInSAR cannot contribute sensible information for sub-centimeter surface displacements. Nonetheless, different zones of sub-centimeter deformation are revealed, but the decrease of accuracy is clearly linked to the vicinity of available GPS measurements by which the joint inversion is constrained.

In general, the north-south component shows a much smoother velocity field than for the other components (Figure 5a). This is because the DInSAR estimates cannot significantly contribute additional information so that the influence of the GPS predominates the most. This is reflected in the deviations between the input GPS and the SISTEM velocity field, which are the smallest in this directional component, with less than 0.5 mm/year for almost all GPS stations. In areas of LOS displacements with no GPS stations, the WLS solution cannot be sufficiently constrained, which becomes visible in the standard error estimates of the WLS solution, with significantly higher values (Figure 5a’). For those areas the distance-decay function weakens the influence of the distant GPS stations that are used for the decomposition. Their estimates, however, do have an influence on the decomposition and might contradict the actual DInSAR data by being part of another kinematic regime. This is especially visible along the coastal part of the Siracusa bay and the horst structure in the hinterland of Ragusa. Since north-south directed motions cannot be captured from the DInSAR datasets alone, the solution fails to disentangle the DInSAR data in the proper way.

The most obvious case is the southern part of the Scicli-Ragusa fault belt, where GPS stations are absent on the western part of the fault line. As a left-lateral strike-slip fault, the expected main direction of horizontal deformation should be in the north-south direction. The nearest GPS station HMDC, which is located about 4 km east of the fault, shows a 0.4 mm/year westwards and a 0.2 mm/year southwards movement with respect to the Noto reference station. The influence of the GPS measurement is therefore lower on the western side than the contradictory DInSAR estimates. In addition, the information west of the fault is weighted more strongly from the DInSAR in the east-west direction due to the lower distance thresholds for the distance-decay function. As a result, the motion could not be well constrained by the DInSAR data, which is not sensitive to the north-south direction and causes large standard errors of the WLS solution. This partly explains the predominant eastwards motion that does not match the expected north-south direction. Also the high values are most likely unrealistic. Indeed, the area features the highest residuals of up to 5 mm/year between the predicted LOS velocities from the SISTEM decomposition compared to the original ones used as input (Figure 6). It cannot be fully excluded that a part of the signal might be influenced by residual atmospheric noise that can be correlated to orographic effects by the local topography. On the other hand, the clear alignment along the fault and the expected active tectonics in this part of the study area allow us to speculate that the method may pick up the respective tectonics in an unprecedented spatial detail, but fails to estimate the right direction and velocity.

The same rationale applies to the Avola fault and its secondary structures in the south-east of the study area. A westward movement of about 2–4 mm/year has been detected which is bounded along the fault. Due to the absence of GPS measurements on the eastern side of the fault, the quantitative and directional displacements might not reflect the actual deformation, but the fault location is obtained in an unprecedented spatial detail.

The distinctive gradient of southwards motion along the Catania-Gela foredeep south of Mount Etna highlights the importance of the presence of GPS measurements in the vicinity for the retrieval of the north-south component. By excluding the IAZZ station south of Catania, which exhibits a southwards directed motion of about −1 mm/year with respect to the NOTO reference station, the spatial pattern of local shortening disappears within the SISTEM result. Other studies, however, have already described this shortening using additional campaign GPS data for this particular zone [88]. Even though the predicted LOS from this station deviates by 2.6 mm/year for the ascending and 2.9 mm/year for the descending track, the information was crucial to derive a realistic velocity field for this part of the study area.

3.4. The Velocity Field in Its Geological Context

Although outside of the main scope of this study, this section seeks to put the obtained velocity field in a geological context and underline the significance of the obtained results. Most recently, structural, seismological, geodetic as well as combined studies strengthen the idea that the Hyblean plateau acts as an independent crustal unit with respect to the kinematic behaviour of the Nubian plate [70,71,88]. The boundaries are located along the Gela-Catania foredeep, the rift zone along the Malta Escarpment and the Scicli-Ragusa line. The retrieved deformation field depicts displacements on all three parts of the plateau as well as various local zones of vertical motions that will be discussed below.

3.4.1. Northern Rim of the Hyblean Plateau

The northern rim of the Hyblean plateau resides in a triangle zone between the collisional belt along the Gela-Catania foredeep, and a N-S oriented rift zone towards the east along the Malta Escarpment (Figure 1). While this area was part of an extensional basin until the Lower-Middle Pleistocene, tectonic inversion set in about 0.85 my ago. This process has been indicated by structural analysis of outcrops within the Scordia-Lentini basin [89]. The area of maximum compression in the SISTEM result is located further north with values of −2.4 mm/year along the Gela-Catania foredeep, while only minor velocity gradients of about −1.1 mm/year are detectable over the horst structures bounding the Scordia-Lentini graben (Figure 5a). This fits with prior geodetic studies based on GPS data solely, that could not detect significant GPS strain over the Scordia-Lentini graben [71]. The authors of this study, however, indicate that during tectonic inversion, faults are reactivated selectively, and argue that active faulting has been displaced northwards. Structural analyses confirm this assumption and locate current associated processes within the Catania plain and its northwards bounding structures [90]. In the Catania plain the SISTEM results are in line with prior GPS based geodetic measurements that report shortening rates of circa 2.4 mm/year [71]. More elevated southwards displacements up to −3.4 mm/year are visible farther north towards the foothills of Mount Etna. This is assumed to be associated with the outward thrusting that originates from the basal decollement of Mount Etna, already reported by past DInSAR studies [91,92,93], and coincides with an anticline visible in the vertical component (Figure 5c).

3.4.2. Scicli-Ragusa Fault Belt

The Scicli-Ragusa fault belt extends from the western side of the horst structures at the northern rim of the plateau, southwards to the coast and continues off-shore up to the Pantelleria-Linosa-Malta rift zone [94]. It is considered to play a key role for the understanding of the kinematics of the central Mediterranean as it divides the eastern part of the plateau from the western part. Indeed, geodetic measurements indicate that the area west of the fault moves in agreement with the entire western Sicilian domain. Consequently, a sinistral strike-slip is present and according to Musumeci et al., 2014 [88] amounts to circa 1.2–1.5 mm/year for the northern part of the fault belt, while the southern branch exhibits lower values around 0.5 mm/year. The SISTEM estimates presented in the previous section are dominated by a strong eastwards motion over the horst around the city of Ragusa, which does not match any other geodetic observations or kinematic models that have been published so far. However, the north-south velocity field, despite being associated with large errors, matches the observations from Musumeci et al., 2014 [88] (Figure 7a). Moreover, the Comiso-Chiaramonte fault belt seems to be associated with accommodated tectonic stress. From this consideration it could be speculated that not the entire Scicli-Ragusa fault belt acts as a transfer zone and that part of the collisional force driven by the Nubian-Eurasian convergence superimposes the local geodynamics. The eastwards driven motion of the SISTEM result could then be partly explained by the NW-SE extension of the Pantelleria-Linosa-Malta rift zone, which according to Serpelloni et al., 2007 [95] amounts to about 1.5 mm/year. Low seismic activity and the absence of GPS measurements in this zone do not allow for an in-depth interpretation. Furthermore, the extension of the observed deformation towards the western side needs to be revealed by additional A-DInSAR data.

3.4.3. Ionian Boundary of the Hyblean Plateau

The eastern part of the Hyblean plateau is influenced by a seismogenic rift zone, which extends NNW-SSE oriented along the entire Ionian coast. Active extension takes place off-shore further to the north with a rate of about 3.4 mm/year [96]. Towards the south, part of this rift zone swings into a SE-NW oriented fault belt, consisting of the Avola fault and the Pozallo-Ispica-Rosolino fault belt, exhibiting extensional rates of about 1.5 mm/year estimated by structural analysis [96]. The SISTEM results do not reveal any of the E-W extensional character on the northern part of the Hyblean plateau, which indicates that most of the seismic potential originating from the rift zone is located off-shore. Farther towards the south, fault-aligned displacements are revealed by both horizontal components. The strongest displacement is notable along the Avola fault line. This is a SW-NE directed normal fault with an expressive scarp of about 290 m and is considered to have caused the foreshock (M${}_L$ = 6.2) of the M${}_L$ = 7.4 earthquake on the 11th of January in 1693 [97]. The horizontal velocities from SISTEM indicate a WNW directed displacement of up to 3 mm/year. This is somewhat contradictory to its extensional character and would indicate that it accommodates part of the main tectonic stress originating from the Nubia-Eurasian convergence, which, in turn, is underlined by a slight uplift of the coastal part up to 1 mm/year. Part of this WNW directed displacement continues along linear fault structures towards the centre of the plateau, which is indeed seismically active. Focal mechanisms of this area indicate strike-slip faulting within this direction [88], and a seismic swarm taking place between October 2011 and July 2012 represents the biggest strain release over the whole plateau over the last 10 years [98]. Interestingly, the horizontal velocity field given by SISTEM indicates local variation over the area of this seismic swarm that took place roughly 1 year after the last acquisition of the A-DInSAR time-series (Figure 7b). This provides further indication that the diffuse displacements captured by the SISTEM decomposition over this area are associated with active tectonic driven deformation.

3.4.4. The Presence of Subsidence

The interpretation of the tectonic deformation over the Hyblean plateau is hampered by various other displacements present in the dataset and partly overlaying the above mentioned tectonically driven processes (Figure 5c). Water pumping is a wide-spread practice in Sicily for agricultural and industrial water use. The extraction of ground water leads to a compression of the aquifer which, in turn, induces the subsidence of the surface. Canova et al., 2012 [25] detected areas of subsidence related to ground water exploitation within the Augusta bay and the area around Villasmundo, about 15 km west of Augusta. They used ERS time-series data for the period of the 1990s and decomposed the DInSAR signals of both tracks by neglecting the north-south directed motion. While the extent of the subsidence areas is comparable, their estimates are about twice as high as the ones retrieved by the SISTEM decomposition. The rates of subsidence depicted within our study exhibit up to −10 mm/year in the Augusta bay and range between −2 to −4 mm/year in the surroundings of the village of Villasmundo. While this might show the consolidation of the soil in time, it could be also related to differences of the DInSAR reference point, which is not indicated in the comparative study.

Another zone of subsidence resides in the south of Catania that has not been reported before. Due to the localised extension it is likely being caused by local water pumping. Values amount up to −8 mm/year, which can cause damages to local infrastructures and industrial buildings located within this zone. Other areas of subsidence are seen towards the Gela-Nappe along the Comiso-Chiaramonte fault belt in the west of the plateau and over the south-eastern part. However, the lack of more detailed information does not allow for a definite interpretation. While tectonic induced stress is speculated along the Comiso-Chiaramonte fault belt indicated by the horizontal velocity field, the revealed subsidence might be driven by the gravitational creeping of the alluvial soils in this zone. In the south-eastern part, the subsidence does also overlap with horizontal displacements, and thus might be partly caused by active tectonic deformation. On the other hand, this area has been reported to be affected by diffuse karst processes [24] that may lead to subsidence, too.

4. Conclusions

In this paper we describe a complete workflow for combining advanced time-series DInSAR with GPS data for the retrieval of a full 3-dimensional velocity field of the Hyblean plateau in south-east Sicily. Due to expected small deformations we applied a data-driven approach and corrected the tropospheric InSAR noise using auxiliary data rather than spatio-temporal filtering. We used the TRAIN toolbox, which is fully compatible with the StaMPS data structure, to estimate a tropospheric correction using the ERA-Interim weather model.

An additional post-processing procedure that enhanced the mean LOS velocities with a gap-filling and consequent multi-directional filtering technique was proposed and proved to be successful. The spatially-limited gap-filling resulted in an increase of coverage and of overlapping pixels which is favorable for further processing of multiple datasets. Due to the multi-directional filter, pixels with increased phase noise were filtered leading to a smoother deformation signal in space, while steep gradients along faults were not affected. Again, this is a straightforward approach that can be easily implemented in post-processing procedures.

For the first time, we applied the SISTEM approach for interseismic deformation retrieval. In accordance to prior studies, the approach seems to work well in areas where a sufficient density of accurate GPS data is available. Predicted LOS velocities from the SISTEM output show generally low residuals of about ±1.5 mm/year and less with respect to the original ones in areas of subsidence and uplift. Despite these errors, patterns of detailed horizontal displacements are revealed and are in line with prior studies for various zones of the plateau. Considerable improvements can be expected by incorporating data from the recently launched SAR satellites, such as Sentinel-1, which allow for much higher temporal sampling and offer a more favorable configuration for A-DInSAR time-series approaches. Together with additional GPS measurements, our approach should be capable of providing more reliable estimates of the actual velocity field and thus allow for an improved understanding of the ongoing kinematics as well as the regional seismic hazard. Given these preconditions, the procedure is transferable also to other tectonic zones. Especially along continental plate boundaries, where the tectonic stress is often distributed over a diffuse network of active faults, the preservation of the DInSAR resolution of our approach allows for the identification of fault locations.