Semi-automated extraction of Deviation Indexes (DI) from satellite Persistent Scatterers time series: tests on sedimentary volcanism and tectonically-induced motions

We develop a methodology based on satellite Persistent Scatterers (PS) time series and aimed to calculate two indexes which are capable to depict the deviation from a deformation model defined a priori. Through a simple mathematical approach, these indexes reproduce the visual process of identification of trend deviations that is usually performed manually by the radar-interpreter, and guide the prioritization of further interpretation for those areas recording significant variations within their motion history. First tests on semi-automated extraction of the Deviation Indexes (DI) from RADARSAT-1 PS data available over Southern Italy allowed the quantification of tectonically-induced land motions which occurred in February 2005 within the town of Naro, and also the clear recognition of the precursors to mud volcano eruptions which occurred in August 2008 in the village of St. Barbara. For these areas, the information level brought by the DI increases and adds onto that of other PS parameters, such as yearly velocity, standard deviation and coherence. Factors exerting influence on the DI are critically tackled within the discussions, together with the analysis of the potentials of these indexes for monitoring and warning activities of geohazards.

Besides the simple use of the average (annual) rates of motion, local scale applications and analyses of single phenomena frequently can benefit from the information stored within the deformation history of each PS (e.g.Tapete et al., 2012), where ground motions are recorded with millimetre precision at each satellite acquisition.In these cases, the full characterization of the deformation behaviour and the quality of a set of radar targets is often insufficient if only based on the yearly velocities, their standard deviation and the coherence of the available targets.Following this principle, some PSIbased studies of geohazards have brought to light the importance of calculating additional indexes which are able to better characterise the information stored within each PS time series (Cigna et al., 2011;European Space Agency, 2012).The extraction of these descriptors can facilitate the recognition of nonlinear components, accelerations, decelerations and -more generally -any kind of deviation from a trend defined a priori.

F. Cigna et al.: Semi-automated extraction of Deviation Indexes (DI) from PS time series
The aim of our paper is to develop the Deviation Indexes (DI) and discuss their suitability to describe further the patterns recorded within the PS time series.Such indexes reproduce -through a simple mathematical approach and reliable semi-automation -the process of identification of trend deviations which is usually performed manually by the radarinterpreters based on their eyesight.The use of the DI provides the operator with a synthetic indicator that quantifies the deviations recorded within the time series, and also orients the prioritization of further steps of PS interpretation and analysis over those zones recording significant variations within their deformation history.
The mathematical background for the calculation of the Deviation Indexes is described in this paper, along with the critical discussion of the results obtained from the tests carried out on the two sites of St. Barbara and Naro, located in Southern Italy, respectively affected by sedimentary volcanism and tectonically-induced structural deformation.

Algorithm development
Most PSI approaches exploit simple linear functions of phase variation through time to model the deformation components and to extract motion series for reflective targets (the PS) identified over the investigated area (Ferretti et al., 2001;Crosetto et al., 2010).The use of more complex functions (e.g.hyperbolic, seasonal) is often more suitable to model geological processes, as they are usually characterised by accelerated/decelerated behaviours or, more in general, nonlinear variations.Examples of nonlinear processes include, for instance: (i) newly developed subsidence features due to groundwater pumping; (ii) precursors to slope failures; (iii) ground collapses due to dissolution; (iv) land compression as effect of changes in loading conditions after building construction or underground excavations; (v) shrink-swell clays; and (vi) tectonic motions.
To address the above remark, recent PSI-based applications employed quadratic, hyperbolic and seasonal functions to better model land and structural deformation, for instance for the cities of Mokpo in Korea (Kim et al., 2010) and the archaeological monuments of Rome in Italy (Tapete et al., 2012).Multi-interferogram stacking approaches, such as the small-baseline technique, are also well suitable to depict nonlinear deformation components (e.g.Berardino et al., 2002).
Although a priori models are chosen during the processing and the extraction of the time series, some signals deviating from the assumed model can be frequently detected (e.g.Cigna et al., 2011;Chen et al., 2012;Notti et al., 2012;Tapete et al., 2012).This may occur as long as these deviations do not compromise the temporal coherence of the radar targets and cause either their exclusion from the subsequent processing steps, or the ascription of their "non-model" components to other phase terms (e.g. the atmospheric phase screen).Accounting for these remarks and limitations, we developed and tested some post-processing tools capable to quickly identify and quantify the deviations recorded within the time series of huge datasets of PS (which may include even thousands of radar targets), and to "dig up" the targets exhibiting a trend which deviates from a certain deformation model.

Input data
The approach experimented in this work starts from the output of a PSI processing, i.e. a dataset of point-wise radar targets (the PS) for which -point by point -the following information are usually available (Fig. 1): -geographic (longitude and latitude) and/or cartographic (easting and northing) coordinates, and height ( -quality parameters of the provided velocity and height estimates, such as velocity and height standard deviations (σ v and σ h , measured in mm yr −1 and m, respectively); -PS coherence, a dimensionless number ranging between 0 and 1, and quantifying the similarity of the time series with the deformation model chosen during the processing, i.e. the goodness-of-fit of the model with the time series under consideration (0, null; 1, perfect fit).

Computation of the Deviation Indexes
Using the input PS dataset, the following steps are then carried out to create the synthetic indexes: The calculation of this index requires the following steps to be performed (Fig. 2): -a simple linear regression is applied to the displacement records (d i ) of the H interval, and a first order function d H (t) that fits the H portion of the time series is extracted in the slopeintercept form: where the slope v H represents the average velocity recorded during H , and intercept d 0 the value of the motion at t 0 .-the variability of the time series during H is estimated in the form of the standard error of the regression (s).This is calculated by combining the differences between the records d i and the displacement values modelled with the linear function, d H (t i ), for all the N H acquisitions included in the H interval: -the function d H (t) is then used to predict the displacement trend of the U interval, by assuming that no trend deviations occured during U with respect to H : -the differences i between the displacement records d i over U and the corresponding predictions at time t i are calculated as: where s is the standard error of the regression, calculated by means of Eq. ( 2).This index is dimensionless and takes on only positive values.The higher is its value, the stronger is the deviation recorded during U with respect to its prediction based on the pattern of H . On the other hand, the closer is DI 1 to 0, the weaker is the trend deviation.The value of DI 1 may be associated to the ratio between the average variation during U and the variation during H , hence a DI 1 of +3.0, for instance, suggests that the average deviations recorded during U are 3 times higher than the variability of H .
b. Approach 2: The second typology of index, DI 2 , is designed for rapid displacements which occurred at certain dates and are expected to have influenced the time series locally, i.e. in the form of a sudden change recorded within the displacement series, such as those occurring during tectonic events or structural collapses.Estimation of this index requires the following steps to be carried out (Fig. 3): -a simple linear regression is applied to the displacement records (d i ) of both the H and U intervals (separately), and two first order functions d H (t) and d U (t) are extracted in the slopeintercept form, respectively: where the slopes v H and v U correspond to the average velocities recorded in the H and U intervals, respectively.Values of v H and v U are not necessarily the same and can differ from each other, especially if an acceleration/deceleration also occurred during U with respect to H . -the Deviation Index DI 2 is then calculated as the difference between the values assumed by the two functions d H (t) and d U (t) at the break time (t b ) of the displacement series: This index is measured in millimetres, and the closer to 0 its values are, the lower is the influence exerted by the occurred event on the displacement time series.On the other hand, significant steps recorded in the time series produce DI 2 of a few to tens of millimetres, corresponding to the displacement recorded at t b .Positive DI 2 indicate steps recorded in the direction towards the sensor, while negative DI 2 are due to rapid motions away from the sensor.It is worth noting that the value of this index is very sensible to the choice of the breaking time t b , and can be misleading if there is no clear trend difference between H and U .

Case studies in the Caltanissetta basin
We tested the above described approaches on the sites of St. Barbara and Naro, in Southern Italy, both located at the edge of the Apennine-Maghrebian thrust belt within the Caltanissetta basin (Fig. 4), a dynamic foredeep basin dating back to a period spanning from the Late Miocene to the Quaternary.The Tortonian clays of the Terravecchia Formation are covered by the evaporates of the Upper Miocene Gessoso-Solfifera Formation consisting of, from the bottom to the top, bedded diatomites of the Tripoli Formation, limestones and gypsum interbedded with clays.The Lower Pliocene marls and marly limestones of the Trubi Formation lie on the Gessoso-Solfifera Formation, are interlayered by the Early Pliocene clays, and covered by Quaternary alluvial deposits (e.g.Lickorish et al., 1999).
For St. Barbara and Naro we used two different sets of PS data obtained from satellite SAR imagery acquired with 24-day revisit time along the sun-synchronous ascending passes of the RADARSAT-1 satellite, from an altitude above the Earth of ∼ 800 km and using the right-looking geometry of acquisition.The RADARSAT-1 radar sensor employs a C-band signal with wavelength λ = 5.6 cm and frequency f = 5.3 GHz, and HH polarization.The employed imagery was acquired in standard beam mode 3 (S3), hence the scenes   have a nominal ground-range resolution of ∼ 30 m, and are acquired using a look angle θ of ∼ 35 • , measured between the LOS of the satellite and the vertical direction.Tilt of the satellite orbit with respect to the N-S direction is equal to ∼ 10 • at the latitudes of the two test sites (i.e.37 • N), hence the LOS is also tilted of the same angle with respect to the W-E direction, i.e. it is approximately WSW-ENE oriented.
The SAR stacks acquired over the two test sites were processed by the Italian company TRE S.r.l. with the Permanent Scatterers InSAR (PSInSAR) technique, by employing the amplitude dispersion criterion for the selection of the targets, and a simple linear model of phase variations through time to extract their deformation components.The PS coherence thresholds determining the acceptance/rejection of the targets during the processing were calculated, for the two stacks, by choosing a false alarm rate of 10 −5 , i.e. assuming that only 1 point every 10 5 could be considered unreliable (Ferretti et al., 2001).The NASA Shuttle Radar Topography Mission (SRTM) Digital Elevation Model, with spatial resolution of approximately 90 m and a linear vertical absolute error of ∼ 6 m (average error for 90 % of the data over Eurasia; Rodríguez et al., 2005), was initially employed during the processing for the subtraction of the topographic phase components.The precision of the topographic information was then improved during the PSInSAR processing itself to get PS elevations as precise as 1-1.5 m for both the sites.For the two analyses, the locations of the reference points (zero deformation points) to which all the motion estimates are spatially related, were chosen within zones assumed as devoid of ground motions, accounting for both their geological and geomorphologic settings, and the coherence distribution of the targets observed over the two stacks during the processing.

Mud volcano eruption in St. Barbara
The first area of interest includes the site of Vulcanelli in the Terrapelata quarter, located S of the village of St. Barbara, in the south-eastern sector of the city of Caltanissetta.The village developed during the XX century after 1930s and, geologically, is built upon the Pliocene Argille Brecciate IV Formation (clays with olistostromic intercalations), which widely outcrop in the urban area and its surroundings (Vallone et al., 2008).
On 11 August 2008 this area was affected by a violent paroxysmal eruption of a mud volcano, which caused damages to urban infrastructure located up to distances of 2 km from the main eruptive vents (Madonia et al., 2011).About 9500 m 3 of grey muddy clays containing hot salty water, methane and lithic fragments of centimetre diameters, were emitted from twenty-one vents due to the high pressure of the gas accumulated below the surface.The expelled muddy deposits reached several tens of metres of height, and covered a circular area of ∼ 12 000 m 2 , with maximum sediment thickness of 3.5 m observed close to the main vents, and minimum thickness of 0.3 m observed at the boundaries of the ejected deposits.Geothermometric analysis carried out by Madonia et al. (2011) located the source of the fluids originating the mud eruption at about 3 km depth.The paroxysmal event lasted several minutes and was anticipated by a telluric event occurred a few hours before in the whole Terrapelata quarter and, contemporaneously, in the neighbouring area of St. Anna.The latter caused opening of fractures on the ground and severe damages to local buildings and hydraulic facilities.
Last documented eruptions that occurred in St. Barbara date back to the late XVIII and early XIX centuries, and are associated to the onshore sedimentary volcanism affecting Sicily.This phenomenon occurs over the accretionary wedge developed in front of the Maghrebian trust belt and originates in the sediments deposited in the Caltanissetta Basin.The geodynamical evolution of the stress field of Central Sicily controls considerably this process.Other areas affected by similar processes include the Maccalube of Aragona, Fuoco di Censo in Bivona, Bissana in Cattolica Eraclea, Marianopoli in Caltanissetta, and Salinelle St. Biagio and Salinelle Stadio in Paternò (Etiope et al., 2002).
To depict the deformation scenario before the paroxysmal event of August 2008 and verify the potential occurrence of precursors, we carried out a satellite study using a dataset of RADARSAT-1 PS derived from 56 acquisitions spanning the interval 13 July 2004-9 August 2008, i.e. till the last available acquisition before the eruption (Fig. 5).The reference point for the PSInSAR dataset was located within the city of Caltanissetta (WGS84 geographic coordinates: 14 • 03 42 E, 37 • 30 00 N; 2 km E of Fig. 5a).
We considered a subset of 315 PS in the surroundings of the mud volcano.Velocity standard deviations σ v for this subset range between 0.38 and 0.99 mm yr with respect to the linear deformation model ranges between 0.66 and 0.97, and the value that occurs more frequently (i.e. the mode or modal value) is centred at 0.70.About 35 % of the targets show coherence higher than 0.80, hence their time series fit well with the linear model.More than 65 % of the PS subset is included within the range of ±2 mm yr −1 , and about 75 % of them show even lower velocities (range of the ±1 mm yr −1 ; Fig. 5a).Most of the targets that show velocities exceeding ±2 mm yr −1 are concentrated in the areas close to the volcano, in the Vulcanelli sector.Maximum velocities are observed to the south (−15.3 mm yr −1 ) and west (+14.3 mm yr −1 ) of the mud volcano, respectively.
We applied the first typology of post-processing analysis and calculated the Deviation Index DI 1 by applying Eq. ( 5).The break time t b was set as December 2007 (i.e. 8 months before the event of August 2008), to consider a sufficiently long period over which a representative set of i could be calculated.The resulting interval H includes 43 scenes, while U consists of the last 13 acquisitions of the data stack.The resulting of DI 1 range between +0.2 and +10.3, with more than 55 % of the targets characterised by DI 1 lower than +1.0 (dark green PS in Fig. 5b), which also coincides with the modal value of the whole observed DI 1 .
A significant amount of PS showed high Deviation Indexes, and their time series revealed the occurrence of strong deviations with the respect to the history of the preceding years.11 % of the PS is characterised by DI 1 higher than +2.0 (orange and red PS in Fig. 5b), indicating that the average deviations recorded during December 2007-August 2008 (i.e.U ) are more than 2 times higher than the variability of July 2004-December 2007.The time series of these PS -such as those of PS A2RJG and A2O8D -showed up to 3-5 cm of progressive LOS movements accumulating just before the event in the direction towards the satellite (Fig. 6).Other time series recorded significant accelerations in the direction away from the satellite, e.g.PS A2NFY.It is easily observable that the targets showing uplift are those located on the western side of the volcano, whilst PS showing motions away from the satellite are located in the southern sector of the volcano.The observed movements likely indicate an expansion of the crater area that preceded the eruption of August 2008.Due to the orientation of the RADARSAT-1 LOS, this expansion is seen as a negative motion in the southern side of the volcano, while it is recorded as a positive motion on the western side of the volcano, where the LOS sees the expansion as a displacement towards the sensor.

Tectonically-induced deformation in Naro
The second case study concerns the urban area of Naro, which was affected on 4 February 2005 by a deformation event likely attributable to tectonics.
The town is built on the south-western flank of a hill made of Upper Pliocene calcarenites and sands dipping to the SW with an average slope of 20 • , and overlaying Middle Pliocene clays (Cigna et al., 2011).The geological setting of this area is controlled by the geodynamic evolution of the Gela Nappe which induces E-W and NW-SE faulting and fracturing in this portion of the Caltanissetta basin (e.g.Lickorish et al., 1999).Indeed, Naro hill shows several well-marked slope changes, morphologic steps and scarps of probable morphotectonic origin oriented along the NW-SE direction, which likely represent the surface evidence of deep-seated tectonic structures.
The event of 2005 was identified within the town with the sudden re-opening of several sub-vertical fractures of centimetre to millimetre aperture.These were concentrated mainly over the buildings and the road infrastructure of the historic old quarter in the north-eastern sector of the town (Vanelle St.) and, secondarily, in the south-eastern quarter.Previously documented occurrences of similar events in the historic quarter date back to the XVII, XIX and XX centuries, and suggest that the recent event of 2005 represented the reactivation of an old phenomenon.On site inspections and field checks carried out in 2005 revealed that most of the fractures were aligned along the NW-SE direction, i.e. parallel to the direction of the major scarps and morphologic steps observable within the urban area (Cigna et al., 2011).Further evidences on the ground suggested that the opening of these fractures was the result of differential sub-vertical motions undergone by the different sectors of the built-up area which are bounded by the above morpho-tectonic discontinuities.To test the potential of our DI, we analysed a dataset of 778 RADARSAT-1 PS covering the monitoring interval March 2003-May 2007, hence including the date of the event of 2005 (Fig. 7).This set corresponds to a spatial subset of a wider PS dataset which includes the whole urban area and its surroundings, and whose interpretation is comprehensively described by Cigna et al. (2011).
Within the selected subset, the velocity standard deviations σ v of the targets range between 0.13 and 0.59 mm yr −1 , and yearly velocities are generally lower than −5.0 mm yr −1 in the central sector of the town, with peaks of −11.0 mm yr −1 within the old historic quarter.PS coherence ranges between 0.57 and 0.97 (with modal value 0.62), and a good portion of the targets fits quite well with the linear model (∼24 % PS shows coherence higher than 0.80).The reference point of the entire processing area was purposely positioned at the location of 13 • 46 55 E, 37 • 17 52 N (WGS84 geographic coordinates), as shown in Fig. 7a.
The second typology of Deviation Index was considered as more appropriate in the case of Naro, since the alteration of the series to be recognised was expected as being a rapid displacement occurred at a certain date.By fixing the break time t b at January 2005 (i.e.just before the event), the H interval resulted to be populated by 23 scenes, while the U by 31 images.Equation ( 8) was then employed to calculate the Deviation Index DI 2 of the PS dataset with respect to the event of 2005.
The resulting DI 2 range between −9.0 mm and +9.3 mm, and their modal value is −1.0 mm.Time series with high DI 2 are those recording a sudden change of deformation behaviour at the time of the event (Fig. 8).This change is quantified by the DI 2 and corresponds to values up to 9.0-9.3mm of displacement occurred between January and March 2005, as a motion either towards or away from the satellite sensor.The spatial distribution of DI 2 shows that different deformation behaviours are recorded for the different blocks bounded by the major NW-SE scarps identified within the  town (Fig. 7b), as already observed by the previous PS study of Cigna et al. (2011).PS moving away from the sensor, such as PS A1JZ8 and A1HN4, are those located in the central portion of the town.On the other hand, the PS moving towards the satellite are concentrated in the SW sector of the analysed area, especially S of the main morpho-tectonic scarp running south of F. Crispi Sq., but also in the northern portion of the town, i.e.NE of Vanelle St (Fig. 7b).
By using the entire PS dataset available for the urban and sub-urban area of Naro, Cigna et al. (2011) performed a manual classification of all the PS time series and discriminated two different classes of targets: affected and unaffected by the event occurred in February 2005.The class of affected PS included targets whose time series were influenced by the event (in terms of either temporary or permanent change in their trend), while unaffected targets included PS whose time series were not influenced by the event (either stable or purely linear time series).Independently from the results of this manual classification, our experimentations with the Deviation Index DI 2 clearly highlight that this difference in the deformation behaviour of the different sectors of the town is already detectable by means of the semi-automated extraction of the DI 2 .Complementarily to the manual and laborious classification performed by Cigna et al. (2011), our DI is able to convert a qualitative classification of the targets (affected vs. unaffected) into a quantitative estimation of the amount of the recorded deviation -expressed in terms of millimetres of motion straddling the event of 2005.Additionally, the sign of DI 2 (i.e., + or −) is an indicator of the direction of the occurred motions.This eases the identification of the sectors of the town which moved towards/away from the sensor.

Comparison of DI 1 and DI 2 with other PS parameters
The Deviation Indexes can effectively facilitate the process of understanding of the temporal and spatial deformation patterns recorded over the observed areas, thanks to the properties of (i) selectivity; (ii) semi-automation; and (iii) reduced time-consumption.All these technical achievements have remarkable impacts on the radar-interpretation of PSI-based deformation evidences.
It is worth discussing in detail the results of the tests carried out on the areas of St. Barbara and Naro, by analyzing the relationships between the values of the two DI, and the three main parameters characterizing the PS time series, i.e. yearly motion velocity, PS coherence, and velocity standard deviations.
i. Deviation Index vs. LOS velocity: the results for St. Barbara and Naro show that, although deviation within the deformation series may be found for targets characterised by high motion velocities, there is not a unilateral relationship between high motion velocities and the occurrence of trend deviations (Figs.9a and 10a).In other words, radar targets showing high annual motion velocities v, not always are those recording significant deviations within their time series.
Indeed, not all the PS moving at high rates in St. Barbara register a deviation (Figs. 5 and 9a).It can be observed that many targets with velocities exceeding ±5 mm yr −1 do not record any deviation during the months preceding the event of August 2008 (compare Fig. 5a and b).
On the other hand, the results for Naro show that strong trend deviations can be found not only for PS moving at high rates (Figs.7 and 10a), but also for targets which apparently are stable (e.g.green PS A1KRY and A1HN4; Figs.7a and 8).
In these cases, the sole value of the average velocity should not allow the identification of the deviations, since the latter are "hidden" by the process of averaging over the entire monitoring interval that -to a first approximation -seems to be stable.This evidence is particularly relevant for the interpretation of PS data and the characterization of the unstable areas.Indeed, it suggests that the radar-interpreter should not overlook those targets which seem to be stable based on the simple observation of their velocities.These targets may have recorded important information which are stored within the time series, but smoothed by the long-term deformation history of the target.
ii. Deviation Index vs. Velocity Standard Deviation: the standard deviation of the PS velocities quantifies the degree of variability/dispersion of the time series, and helps to define the precision of the velocity estimates.
The lower the standard deviation, the higher is the precision of estimated velocity.PS standard deviation increases moving away from the reference point, and generally decreases when more populated and longer image stacks are used during the PSI processing.Evidences similar to those noticed for the PS velocities, are found by comparing the two DI with the velocity standard deviations.
For St. Barbara and Naro, it is found that, although high DI usually correspond to high standard deviations σ v , not all the PS with high velocity standard deviations return high values of deviations (Fig. 9b and 10b).
Again, the level of information brought by the DI adds onto that of the velocity standard deviations, and helps to discriminate PS with high variability of their time series from those with clearly distinguishable trend deviations.
iii.Deviation Index vs. Coherence: PS coherence quantifies the degree of correspondence of the time series with the model assumed during the processing, and gives an idea of the inconsistency of the motion history with respect to the employed deformation model.It is clear that any deviation from the expected trend will have a significant role on the value of coherence and will cause it to decrease.Interferometric phase noise will negatively affect the values of PS coherence, as well.
Cross-comparison of the estimated DI with the PS coherence for St. Barbara and Naro (Figs. 9c and 10c) highlights that high values of DI are found for PS with low coherence.However, the fact of having a low coherence does not always means for a target to be

Factors influencing the DI
The parameter exerting the major influence on the DI experimented in these tests is the breaking date (t b ).Its choice plays a fundamental role during the estimation of the Deviation Indexes, because it (i) controls the length of the intervals H and U , and the reliability of their recalculated velocities; and (ii) acts on the possibility of capturing (or not) the occurred deviations that might have happened at any time during the monitored interval, and thereby being detectable (or not) as a result of the selection of a certain t b .The breaking date should consequently be calibrated by accounting for the typology and characteristics of the deformation expected in the monitored area and their hypothesised behaviour, by choosing a date that reasonably represents a potential time for the area to undergo a deformation change.
A priori factors that influence the DI include the model employed during the PSI processing to unwrap the deformation phases, and the interferometric aliasing effects.The deformation model may smooth a priori the time series, and may induce any signal that highly differs from the model to be attributed to other phase components, such as noise or the atmospheric phase screen.Further risk can arise for those PS showing a history of deformation that differs significantly from the model, since they might be given very low coherence or even be discarded during the processing itself.
As mentioned above, influence from aliasing effects is not negligible as well.When motions between two successive acquisitions exceed a quarter of the signal wavelength, λ/4 (e.g. 14 mm for C-band sensors, such as that onboard RADARSAT-1), aliasing effects and phase-unwrapping errors may occur (Hanssen, 2005).Smoothed time series can be the result of phase unwrapping errors, and the motion recorded at a certain date is seen by the technique as lower than it actually is.This thereby causes underestimation of motion deviations, and consequently the magnitudes of the detectable deviations are lower than the actual ones.

Operational exploitability and perspectives
The use of DI contributes towards the reduction of time-consumption with regard to manually-performed time series analysis.Nevertheless, as previously discussed, (semi-)automated computations can lead to an occasional loss of information, if the typology and the input parameters of the employed DI are not sufficiently suitable to describe the phenomena affecting the observed areas.Manual checks of the time series of PS subsamples might be a sustainable validation method to assess the reliability of the automated computations, under the input parameters (e.g.type of DI, choice of t b , length of H and U ) and in face of the phenomenology of the observed natural and/or human-induced instability processes.
Further elements to be taken into account are the specific objective of the analysis and the scale at which we want to perform it.With regard to damage assessment over buildings and single structures due to, for instance, volcanic and tectonic events, the relevance of each extracted DI is to be ascertained after the corresponding PS are verified as being actually related to the object(s) of interest.Relying only on the results of semi-automated extraction of DI might increase further the risks of PS misinterpretation, as already discussed for the so-called "unrelated PS" by Tapete and Cigna (2012).
The simplicity of the mathematical principles behind the Deviation Indexes DI 1 and DI 2 facilitates their export to other case studies.The calculation of the DI can be consequently performed by a wide spectrum of radar-interpreters, also including non-expert users, and can become a routine post-processing step.
Looking at potential future uses of the DI, the improved temporal sampling of PSI data that can be achieved by using new and upcoming SAR missions (e.g.TerraSAR-X, COSMO-SkyMed, or the Sentinel-1 constellation) encourages further development of these semi-automated approaches for the advanced exploitation of PS time series.Implementation of indexes such as DI 1 and DI 2 needs to account for the peculiarities of new high resolution X-band satellite data, and the high volumes of the derived products of the upcoming C-band Sentinel-1 constellation.In this regard, Bovenga et al. (2012) recently found that the densities of the PSI datasets retrieved through the processing of SAR imagery at high resolution (COSMO-SkyMed StripMap HIM-AGE) can be up to one order of magnitude higher than that achievable by means of medium resolution data (EN-VISAT ASAR Image Mode IS2).By processing TerraSAR-X StripMap images, Crosetto et al. (2010) even retrieved PS densities 47 times higher than previous PSI processing of medium resolution imagery.Bearing these perspectives in mind, semi-automated extraction of DI offers a practical response to face the availability of larger volumes of PS datasets and related time-consumption issues.This advantage is also essential in relation to increasing target densities provided by new processing algorithms (e.g. the SqueeSAR as second generation of the PSInSAR technique), as confirmed by several authors (e.g.Notti et al., 2012;Tapete et al., 2012).
The weekly temporal samplings of the new satellite missions suggest that the time has come for the scientific and operational communities to design and develop reliable and robust approaches for the semi-automated analysis of PS time series.As recently discussed by Mazzanti et al. (2012), any advanced tools for the analysis of satellite series can support monitoring and warning activities in the framework of hazard and risk management.The future development of (semi-)empirical forecasting methods which employ long satellite series of ground motions for the identification of hazard (re)activations will likely have a remarkable role during this decade.Such typology of applications may involve, for instance, the prediction of landslide failures for phenomena characterised by time-dependant strain patterns (Mazzanti et al., 2012), or the setting up of early/early-stage warning (Tapete et al., 2012) systems for the preemptive recognition of precursors to structural collapses of buildings and infrastructure.To this aim, the use of approaches which build upon those for the calculation of the Deviation Indexes DI proposed here will surely contribute to the implementation of such systems and models.

Conclusions
Two different Deviation Indexes (DI) -DI 1 and DI 2 -are proposed to improve current methods of radar-interpretation of PS time series.These indexes, distinguished with regard to the typology of natural hazards at hand, are capable of recognizing deviations of the LOS displacements from the deformation model defined a priori during the PSI processing.
The tests carried out on the case studies of St. Barbara and Naro demonstrated the suitability and effectiveness of the DI to analyse the superficial evidences of ground and structural instability, such as those due to mud volcanoes eruptions and tectonic events.Hence the use of the DI can be applied for event-driven impact assessment.The results over St. Barbara revealed DI 1 higher than +2.0 for 11 % of the analysed PS, thereby indicating that the average deviations recorded during the eight months preceding the event of August 2008 were more than 2 times higher than the variability of the times series over the H period July 2004-December 2007.Up to 3-5 cm of progressive LOS movements accumulating just before the event were observed for the time series with higher DI 1 .
The benefits obtained for the case of Naro consisted of the identification of the buildings actually affected by the event of 2005, associated with the quantification of the displacements occurred between January and March 2005, i.e. the period straddling the event of February.The spatial distribution of the higher DI 2 (up to 9.0-9.3mm) also highlights a good level of agreement with the different blocks that are bounded by the major NW-SE scarps identified within the town.This zoning based on the semi-automated extraction of the DI 2 confirmed the corresponding manual analysis published by Cigna et al. (2011), with a significant improvement in terms of time required for the PS classification and trend deviation detection.The comparison with the previous analysis also provided validation of the semi-automated method proposed here.
As being complementary to the other parameters commonly used to describe PS datasets, DI 1 and DI 2 are expected to support PSI-based deformation analyses for natural hazard management.Impacts of routine implementation of the DI might also be perceived during emergencies, when the rapidity of the automation is essential to retrieve an immediate indication of the susceptible/affected areas, and to address further analysis of the data and plan targeted field surveys.In these terms, the test performed over St. Barbara is of particular relevance, since sedimentary volcanism and seismicinduced mud volcanoes are not so rare and can have severe impacts on buildings and infrastructure, as occurred during the earthquake of May 2012 in Emilia Romagna region, Italy.
The added value of these DI lies in their practicality to facilitate the trend deviation identification while handling huge volumes of PS data (up to several thousands of targets per square kilometre), as those being increasingly provided after PSI processing of SAR imagery from new high resolution satellite sensors and/or latest multi-interferogram processing algorithms.

Fig. 1 .
Fig. 1.Methodology employed for the identification of trend deviations within PS time series.

Figure 3 .Fig. 3 .
Figure 3. Estimation of the Deviation Index DI 2 through the assessment of the step record 621 between the intervals H and U. 622 623 24

Figure 4 .
Figure 4. Location of St. Barbara and Naro over the schematic structural setting of Sicily (Italy).

Fig. 4 .
Fig. 4. Location of St. Barbara and Naro over the schematic structural setting of Sicily (Italy).

Figure 7 .
Figure 7. Naro test site: LOS velocity (a) and Deviation Index DI2 (b).The reference point is 39 represented as a green star.40

Fig. 7 .
Fig. 7. Naro test site: LOS velocity (a) and Deviation Index DI 2 (b).The reference point is represented as a green star.

re 9 .Fig. 9 .
Fig. 9. Comparison of Deviation Index DI 1 with LOS velocity (a), its standard deviation (b) and PS coherence (c), for the PS time series of the test area of St. Barbara.

Figure 10 .
Figure 10.Comparison of Deviation Index DI 2 with LOS velocity (a), its standard deviation (b) 657 PS coherence (c), for the PS time series of the test area of Naro.658 Fig. 10.Comparison of Deviation Index DI 2 with LOS velocity (a), its standard deviation (b) and PS coherence (c), for the PS time series of the test area of Naro.

F.
Cigna et al.: Semi-automated extraction of Deviation Indexes (DI) from PS time series Semi-automated extraction of the DI from RADARSAT-1 (2003-2008) PS time series allowed the recognition of the precursors to the eruption occurred in St. Barbara in August 2008, and the quantification of land motions occurred in February 2005 within the town of Naro.
LOS) estimates of annual motion velocity (v, measured in mm yr −1 ), corresponding to the average rates of deformation observed over the entire monitored interval [t 0 − t n ];n records of LOS estimates of ground displacement (d i , measured in mm) occurred during the monitored interval [t 0 − t n ], and recorded acquisition by acquisition; h, measured in m); Nonlin.Processes Geophys., 19, 643-655, 2012 www.nonlin-processes-geophys.net/19/643/2012/ -line-of-sight ( Estimation of the Deviation Index DI 1 through the comparison of pattern U with its prediction derived by using pattern H . 1.A breaking time, t b , is arbitrarily chosen within the monitoring interval [t 0 −t n ].This date represents the position at which the time series of the available PS will be cut, and allows the distinction of two different subintervals within the time series, i.e. the historical (H ) interval [t 0 − t b ] and the updated (U ) interval [t b − t n ].These are the two time ranges that will be compared to identify any potential deviations occurred during the monitoring period.2.The monitored period [t 0 −t n ] is then sub-sampled using the break time t b , and the two different sub-intervals H and U are consequently identified within each time se-