Investigating Spatial Patterns of Persistent Scatterer Interferometry Point Targets and Landslide Occurrences in the Arno River Basin

Persistent Scatterer Interferometry (PSI) has been widely used for landslide studies in recent years. This paper investigated the spatial patterns of PSI point targets and landslide occurrences in the Arno River basin in Central Italy. The main purpose is to analyze whether spatial patterns of Persistent Scatterers (PS) can be recognized as indicators of landslide occurrences throughout the whole basin. The bivariate K-function was employed to assess spatial relationships between PS and landslides. The PSI point targets were acquired from almost 4 years (from March 2003 to January 2007) of RADARSAT-1 images. The landslide inventory was collected from 15 years (from 1992–2007) of surveying and mapping data, mainly including remote sensing data, topographic maps and field investigations. The proposed approach is able to assess spatial patterns between a variety of PS and landslides, in particular, to understand if PSI point targets are spatially clustered (spatial attraction) or randomly distributed (spatial independency) on various types of landslides across the basin. Additionally, the degree and scale distances of PS clustering on a variety of landslides can be characterized. The results rejected the null hypothesis that PSI point targets appear to cluster similarly on four types of landslides (slides, flows, falls and creeps) in the Arno River basin. Significant influence OPEN ACCESS Remote Sens. 2014, 6 6818 of PS velocities and acquisition orbits can be noticed on detecting landslides with different states of activities. Despite that the assessment may be influenced by the quality of landslide inventory and Synthetic Aperture Radar (SAR) images, the proposed approach is expected to provide guidelines for studies trying to detect and investigate landslide occurrences at a regional scale through spatial statistical analysis of PS, for which an advanced understanding of the impact of scale distances on landslide clustering is fundamentally needed.


Introduction
Landslides are a major type of natural hazards around the world.They may pose a great threat not only to human beings but also to infrastructure objects, leading to enormous losses in lives and economics.Italy is a country very susceptible to landslide hazards.For example, during the 20th century, landslides caused 5831 deaths and the homelessness of more than 700,000 people in Italy [1,2].In addition, each year in Italy, landslides can generate direct economic losses of 1-2 billion Euros, representing 0.15% of the national gross domestic production (GDP) [3].
Persistent Scatterer Interferometry (PSI) is an Interferometric Synthetic Aperture Radar (InSAR) approach processing multiple interferograms from temporal Synthetic Aperture Radar (SAR) images.The outputs of PSI, namely PSI point targets, often simplified as Persistent Scatterers (PS), are point-wised benchmarks having long-term stability in phases.PSI point targets can be achieved from various approaches of SAR images processing.For instance, PSInSAR TM technique was developed by Ferretti et al. [4,5] and further improved by Colesanti et al. [6] to assess phase stability of the scattering barycenter with low differential atmospheric contribution regardless of normal and temporal baseline.SBAS (Small Baseline Subset), presented by Berardino et al. [7], Lanari et al. [8] and Casu et al. [9], is based on interferograms of small baseline and uses the Singular Value Decomposition (SVD) approach to combine unwrapped differential InSAR interferograms.IPTA (Interferometric Point Target Analysis) introduced by Werner et al. [10] and Strozzi et al. [11] is able to obtain stable benchmarks from low coherent regions with phase estimation from large baseline.CPT (Coherence Pixel Technique), as proposed by Mora et al. [12] and Blanco-Sanchez et al. [13], can estimate both linear and nonlinear components of ground deformation.StaMPS (Stanford Method for Persistent Scatterers), as reported by Hooper et al. [14,15], particularly employs the spatial correlation of interferogram phases to detect persistent radar benchmarks.STUN (Spatio-Temporal Unwrapping Network), as recommended by Kampes [16], enables a combination of displacement models with spatial networks for phase unwrapping through single-master interferograms.SPN (Stable Point Network), as designed by Crosetto et al. [17] and Herrera et al. [18], pays particular attention to pixels with interferometric coherence, amplitude stability, and spectral consistency for PS extraction.In addition, TCPInSAR (Temporarily Coherent Point InSAR), as suggested by Zhang et al. [19], involves temporarily coherent point networking and a least squares estimator, and is capable of extracting deformation signals without the need for phase unwrapping.
In recent years, PSI has increasingly shown its usefulness in landslide studies in various ways.Firstly, several case studies have utilized PS for detection and mapping of landslides.For example, Bovenga et al. [20] employed PS to detect slope instability in urban and peri-urban areas; Cascini et al. [21] exploited PS to validate and update the landslide inventory at 1:25,000 scale; Righini et al. [22] presented a systemized approach to prepare the landslide inventory from PS; Del Ventisette et al. [23] proposed to use PS information to update regional inventories of slow-moving landslides and deep-seated slope deformations.Secondly, PSI is essential for landslide monitoring which can be fulfilled by combining ground data such as leveling and GPS measurements.For instance, Herrera et al. [18,24] used both C-band ERS and ENVISAT data and X-band TerraSAR-X images to monitor landslides in Spain; Greif and Vlcko [25] utilized PS processed from ERS and ENVISAT images to monitor the post-failure deformation in Slovakia; Tofani et al. [26] combined PS and geotechnical measurements to monitor landslide in Northern Apennines in Italy.Thirdly, PSI is effective for investigating landslide mechanism.For example, Hilley et al. [27] performed one of the earliest studies to analyze landslide dynamics from PS and precipitation data; Colesanti and Wasowski [28] integrated field investigation with acquisition geometry to understand landslide mechanism in Liechtenstein Alps; Calo et al. [29] exploited PS combined with ground measurements to assess the activity of landslides; Zhao et al. [30] attempted to understand mechanism of slope failure and define corresponding rainfall thresholds for landslide early warning purposes by using PS covering California and Oregon, USA.Other representative applications of PSI in landslide studies include: Cascini et al. [31] suggested to project the line-of-sight (LOS) deformation vector obtained from PS to the steepest slope direction supposing that the targeted mass movement is translational landslide; Cigna et al. [32] suggested a framework for assessing landslide intensity from PSI-based matrix approach; Lu et al. [33] proposed a method for quantitatively assessing landslide hazard and risk at regional scale from PS; Bianchini et al. [34] showed how to use PS to improve maps of damage assessement due to landslide occurrence.
Spatial clustering delineates the excess of events or values in geographic space determined by a null hypothesis of the expected spatial pattern, distinguishing geographic events comparing higher or lower spatial densities amongst neighboring areas so as to better judge the data structure and reject the hypothesis of randomness among events [35][36][37].Landslides tend to happen in or nearby places where they have occurred before [38].Landslides appear to be spatially clustered considering that they often occur in places where mass movements were present in the past [39,40].This clustering behavior of landslides was noticed by several literatures focusing on different areas of the world.For example, Trauth et al. [41,42] revealed the effect of spatial clustering of landslides which may document Quaternary climate changes in northwestern Argentina; Jarman [43] indicated that failures of rock slopes tend to cluster in areas in the Highlands of Scotland; Lacoste et al. [44] observed the existence of spatial clustering of multi-scale landslides along coastal areas in New Zealand; Tonini et al. [37] discovered a clustering spatial pattern of landslides in the Rhone Valley of Switzerland.
This spatial pattern of clustering found on landslides was noticed by some studies utilizing PSI point targets for landslide detection and investigation purposes.For instance, Hererra et al. [18] observed that PSI point targets tend to cluster on landslide scarps; Bianchini et al. [45] attempted to detect landslide hotspots where higher hazard level was revealed by high density PS with higher velocity than neighboring benchmarks; Lu et al. [46] reported the clustering behavior of extremely slow-moving landslides in the Arno River basin in Central Italy and proposed a hotspot and clustering analysis of PSI point targets to automatically identify areas with a large number of landslide occurrences; Bianchini et al. [47] tried to develop a systemized approach to evaluate the state of activity of landslides from PS clusters.However, none of these studies has noticed the influence of scale distance in clustering analysis of PS.For example, none of these studies has indicated from which and to which scale distances these clustering analyses of PS can be applied and at which scale distances spatial patterns of clustering, repelling and randomness are significant.Additionally, these studies failed to differentiate the impact of landslide types, PS velocities and acquisition geometry on the suitability of PS clustering analysis for landslide investigation.This calls for a statistical approach which can analyze spatial relationships between different types of landslide occurrences and PS of various velocities and orbits, with particular consideration of effective scale distances.Tonini et al. [37] introduced the univariate K-function to test if landslides mapped in the Swiss Rhone Valley are clustered or randomly distributed.A similar approach was also reported by Tonini et al. [48] which describe the spatial pattern of rock falls in Catalonia, Spain.The univariate K-function has the advantage to characterize spatial patterns at corresponding scale distances, however, it has limitation in analyzing spatial patterns between two different groups of data within the study area; for example, it cannot evaluate the spatial relationships between landslide occurrences and PSI point targets.
To overcome this limitation, the bivariate K-function which is able to describe the spatial pattern between two types of datasets, was deployed to investigate spatial relationships between different types of landslide occurrences and PS with various velocities and acquisition orbits.The Arno River basin in Central Italy was chosen as the study area because crosses the Apennines chain and is very susceptible to landslide hazards, and additionally because of the availability of a complete landslide inventory mapped over 15 years, which documents and classifies different types of mass movements across the basin.The primary purpose is to evaluate whether PS processed from 4 years of RADARSAT-1 images are spatially clustered (spatial attraction) or randomly scattered (spatial independency) on different types of landslides, and to evaluate their scale distances of corresponding spatial patterns.The outcomes are expected to provide a priori understanding of spatial pattern and their corresponding scale distances before the spatial clustering analysis of PS for landslide detection and investigation purposes, so as to help understand which types of landslides can be detected by the clustering analysis of PS and at which scale distances the analysis is effective.In particular, four types of landslides (slides, falls, flows and creeps) classified according to the Hutchinson [39] classification and four categories of PS corresponding to different velocities and orbits (slow-moving ascending, slow-moving descending, extremely slow-moving ascending and extremely slow-moving descending) were concentrated for the bivariate K-function analysis in this study.
The remaining part of the paper is organized as follows: Section 2 presents the study area of the Arno River basin; Section 3 introduces the employed datasets (PS and landslide inventory) for the spatial pattern analysis; Section 4 introduces the employed bivariate K-function analysis; Section 5 demonstrates the derived results for the spatial pattern analysis between PS and landslide occurrences; Section 6 renders discussion on materials, methodology and results; Section 7 concludes the whole study.

Geographic Location
The Arno River basin (Figure 1) is situated in central Italy with 98.4% and 1.6% of the territory in Tuscany and Umbria, respectively.The total area of the basin is about 9130 km 2 , with the hydrographic area of about 8200 km 2 .In the Arno River basin, about 78% (7190 km 2 ) of the territory is situated in mountainous and hilly regions.The mean elevation of the whole basin is about 235 m above the sea level, with 55% of the basin lower than 300 m, 30% between 300 m and 600 m, 10% between 600 m and 900 m, and another 5% higher than 900 m.The basin has about 2.2 million inhabitants which are distributed within 166 municipalities.The major cities include Arezzo, Pisa and Firenze as marked in Figure 1.

Geological Settings
The Arno river basin is situated across the Northern Apennines orogenic belt which is a complex thrust-belt system formed by the overlying of several tectonics units, built during the Tertiary under a compressive regime followed by extensional tectonics from the Late Miocene [49].The latter phase brought several horst and graben structures with an alignment of NW-SE which induced the emplacement of Neogene sedimentary basins that could arise during the Pliocene and Pleistocene [50].
The basin sediments are of marine (west) or fluvio-lacustrine origin (east), depending on the location of the intra-mountainous basins, and are formed from marine transgression and regression cycles (west to east), which occurred during Miocene and Pliocene.Marine transgression and regression cycles occurred during the Miocene and Pliocene while fluvio-lacustrine sedimentation occurred during the Pliocene-Quaternary ages [51].
The drainage system of the Arno River is accordingly influenced by this structure, consequently bringing the prevailing NW-SE trending streams.Four main ridges can be identified:

Landslide Hazard within the Basin
The Arno River basin is very susceptible to hazards such as landslides, floods and subsidence [52].In terms of landslide hazard, the total affected area in the basin is about 800 km 2 .Some 27,503 landslides were mapped from 1992 to 2007 at the scale 1:10,000.The activities of landslides are classified as stabilized, dormant and active, according to the classification proposed by Cruden and Varnes [53].Major factors controlling or changing the landslide activities in the Arno River basin include intensive rainfall, snowmelt, soil erosion and anthropogenic activities, particularly the agricultural mechanization for wine and olive-oil production [54].The detailed information regarding the landslide inventory will be introduced in the next section.Theses landslides have produced great risks to vulnerable elements with the high density population in the Arno River basin: more than 16,000 civil buildings, 460 km 2 of industrial area and 350 km of roads have been affected by landslides [46].The economic losses resulting from landslides in the next 30 years have been estimated to be between 3.22 and 6 billion Euros throughout the whole basin [33,49].

Landslide Inventory
In this study, the landslide inventory, which was initially mapped by the Arno River Authority and subsequently updated by the SLAM (Service for Landslide Monitoring) project of European Space Agency (ESA), was made at the scale of 1:10,000 throughout the whole basin.Due to the large spatial extent, the inventory was primarily created from 15 years (spanning from 1992 to 2007) of remote sensing data and topographic maps, including aerial photos, optical satellite images, regional topographic maps, Digital Elevation Models (DEM), radar images and their derivative PS.The detailed information about the employed datasets is summarized in Table 1.Besides, intensive field surveys were performed as auxiliary approaches for investigating hazardous areas, particularly for validation of mapping results, confirming landslide typology and the state of activity, and identifying unmapped landslides.The mapped landslides were further classified with reference to the classification suggested by Hutchinson [39], chiefly focusing on morphological and hydrological features of mass movements related to geology and hydrogeology for understanding landslide mechanisms.In total, 27,503 landslides were mapped in the Arno River basin.As shown in Table 2, 20,578 of them are rotational and translational slides, accounting for 74.8% of the total mass movements mapped in the inventory.The second most frequently occurring mass movements in the basin are creeps, representing 19.4% of the total number of landslides.The remaining types of landslides are flows and falls, with a total number of 1258 (4.6%) and 330 (1.2%) events, respectively.For each type of mass movement, in order to differentiate the state of activity, the mapped landslides were classified into three categories: stabilized, dormant and active, with reference to the classification proposed by Cruden and Varnes [53] and modified by Farina et al. [54].The mapped landslides were then converted to points datasets by calculating centroids of all inventory polygons for the point-based bivariate K-function analysis and to reduce the computation time under the R environment.

PSI Point Targets
The PSI point targets employed in this study were processed from about 4 years (from March 2003 to January 2007) of RADARSAT-1 images by Tele-Rilevamento Europa (TRE) on behalf of the Arno River Basin Authority based on PSInSAR TM technique [4][5][6].In total, as summarized in Table 3, 102 RADARSAT-1 images, including 54 ascending and 48 descending scenes, were employed for extracting those PSI point targets.These images were acquired with the S3 beam mode of RADARSAT-1 satellite, providing an incidence angle ranging between 30 and 37 degrees.More than 700,000 PSI point targets were derived with a defined coherence level above 0.6.The precision of the estimated displacement rate along line-of-sight (LOS) varies from 0.1 to 2 mm/year, depending on coherences and distances to the reference point.The accuracy of geocoding is within 10 m in range direction and within 5 m in azimuth direction.The PS located in the flat areas (slope < 3°) were then masked so as to only focus on the PSI point targets located in the hilly and mountainous areas.Cruden and Varnes [53] suggested a velocity threshold of 16 mm/year to differentiate extremely slow-moving landslides from very slow-moving landslides.This threshold may be lower for displacement rates estimated by PSI point targets for landslide investigation.This is because the velocity scale proposed by Cruden and Varnes [53] is based on the peak velocity of the mass movement while displacement rates acquired from PS is on the basis of average velocities during the time span of processed SAR images, which may be significantly lower than the maximum velocity [32].In this study, the velocity proposed by Bianchini et al. [45] and Righini et al. [22]-for discriminating between extremely and very slow-moving landslides-for PS was used, which decreases this threshold to 10 mm/year.In addition, considering that PS velocities ranging from 0 to 2 mm/year usually contain noises and the precision of displacement is up to 2 mm/year, the PSI point target identified for very slow-moving events was defined from 4 to10 mm/year in this study.In consequence, with additional consideration of the acquisition orbits of the RADARSAT-1 satellite, four types of PS were focused on for the further bivariate K-function analysis: extremely slow-moving ascending PS, very slow-moving ascending PS, extremely slow-moving descending PS and very slow-moving descending PS.

Bivariate K-Function Analysis
The spatial relationship between PS and landslide occurrences was analyzed through the spatial statistics of the bivariate K-function analysis.Similar to the univariate K-function, the bivariate K-function is on the basis of Ripley's K-function, the object of which is to detect deviations from spatial homogeneity incorporating a similar number of points distributed in a given circular area, particularly judging whether spatial patterns of point datasets tend to be clustered, dispersed or randomly scattered [55,56].Spatial clustering is such deviation that can be found lacking spatial homogeneity through Ripley's K-function using the probability model of the Poisson point process inside a plane with constant intensity function.Nevertheless, compared to the univariate K-function, the bivariate K-function further extends the use of Ripley's K-function into multiple point datasets, capable of evaluating spatial patterns between various groups of point datasets [57].In this study, the bivariate K-function was utilized to investigate the spatial dependence between PS processed from about 4 years (from March 2003 to January 2007) of RADARSAT-1 images and different types of landslides mapped within the Arno River basin.
As a tool for analyzing spatial patterns of point datasets, the bivariate K-function is judged as the expected number of category 1 within a given distance h of an arbitrary point of category 2 [58]: where in this study the subscript i denotes a type of landslide occurrence such as slide, flow, fall or creep, and the subscript j represents a group of PS classified by the acquisition orbit and velocity scale proposed by Righini et al. [22] for landslide classification.N is the total number of the analyzed point datasets in the study area.λ and λ represent the intensity (number per unit area, λ / , A is the total area) of landslide occurrences and PSI point targets, respectively.S i and S j are location of landslide and PS, respectively.The weight function w s , s contributes as the edge correction.
The indicator function s s h judges whether a PSI point target S j is within a distance h of a landslide S i .Furthermore, in order to assist graphical understanding of the outputs of the bivariate K-function analysis, the L transformation [59,60] was further performed, calculated as: where h denotes the outputs of the bivariate K-function within the specific distance h.h characterizes the outcomes of the L transformation on h , within the corresponding searching distance h.It can stabilize the variance and under Complete Spatial Randomness (CSR) the expectation for h is zero [61].Moreover, 99 Monte Carlo simulations were accomplished using h calculated for each simulation, so as to build tests for CSR for judging whether the observed events are constant with a homogeneous Poisson point process.If the output of h is detected outside the envelope, the null hypothesis can accordingly be rejected, suggesting an existing pattern of spatial dependence, for instance, clustering or repelling, between PSI point targets and landslide occurrences.Alternatively, the random spatial pattern between these two datasets can be determined if the output of h is measured inside the envelope.During the testing, the significance level (α) used for rejecting the null hypothesis at each distance h is calculated as: where N s is the number of simulations.As a result, for 99 Monte Carlo simulations, the confidence level of 98% is determined.
The bivariate K-function is thought to be a cumulative function that accumulates spatial information on the distribution of points up to each scale distance.As a result, the bivariate K-function is a distance-based method that estimates average number of neighboring points in a specified region inside a given radius, namely the scale distance, over the density of the whole study area.
h increases along different radius (h), or scale distances, if the average number of neighboring points is higher than the density of the whole study area; on the other hand, h decreases along different radius (h), or scale distances, if the average number of neighboring points is lower than the density of the whole study area [62].In this study, the bivariate K-function counts the number of neighboring PSI point targets for each landslide within a given range of multiple distances from 0 to 80 km. Figure 2 shows the basic concept of this analysis based on multiple distances.In a geographic information system (GIS), this type of analysis based on multiple distances functions practically as operations of "zoom in" and "zoom out".Therefore, alternatively speaking, the bivariate K-function enables a calculation at multiple distances and it reveals how spatial pattern of various point distribution varies with scale.The L transformation of the bivariate K-function is then able to look at clustering from each point with every other point across all distances.This justifies why the distance h is also recognized as "scale distance" in this study.Based on the bivariate K-function, investigation was focused on the spatial patterns between various types of PS related to landslide velocity scale and acquisition orbits (very slow-moving ascending, extremely slow-moving ascending, very slow-moving descending and extremely slow-moving descending) and different categories of landslides mapped in the Arno River basin (slides, flows, falls and creeps).The statistical analysis and graphical visualization were accomplished under the free R programming environment, particularly with the extended packages of Spatsat [63] and Splancs [64,65] for the bivariate K-function analysis.

Slides
The bivariate K-function analysis on spatial dependence between PS and slide occurrences rejected the null hypothesis that PSI point targets appear to cluster similarly across all types of slides in the Arno River basin. Figure 3 shows the result of the bivariate K-function analysis on PS with reference to slides of different states of activities, and the corresponding scale distances of their spatial patterns are summarized in Table 4 with the 98% confidence interval.The result indicated that no significant spatial dependence can be observed between ascending PS and occurrences of active, dormant and stabilized slides in the Arno River basin.In terms of descending PS, spatial randomness was detected between extremely slow-moving PS and occurrences of active, dormant and stabilized slides.However, spatial clustering can be noticed for very slow-moving PS on active and dormant slides, with corresponding clustering distances of 57-58.4km and 0-8.8 km, respectively.No evident spatial clustering or repelling was found for very slow-moving descending PS on stabilized slides.

Flows
The bivariate K-function analysis on spatial relationships between PS and flow occurrences is able to reject the null hypothesis that PSI point targets tend to cluster in similar spatial pattern across all types of flows within the Arno River basin. Figure 4 illustrates the result of the bivariate K-function analysis on spatial dependence between PS and occurrences of active, dormant and stabilized flows, with corresponding scale distances of spatial patterns summarized in Table 4. Similar to the spatial relationship with slides, the ascending PSI point targets exhibited only spatial randomness on active, dormant and stabilized flows.Also, no evident spatial dependence can be observed for extremely slow-moving descending PS across all types of flows in the Arno River basin.Nevertheless, clusterings of very slow-moving descending PS were found significant on active, dormant and stabilized flows, with corresponding clustering ranges of 1.7-2.7 km, 2.1-4.6 km and 73.8-78.1 km, respectively.

Falls
The result of the bivariate K-function analysis on PS and fall occurrences rejected the null hypothesis that PSI point targets tend to cluster in similar spatial pattern across all types of falls in the Arno River basin.The result is visualized in Figure 5, with corresponding distances of spatial patterns summarized in Table 4.The analysis revealed that both ascending and descending PSI point targets, with very slow-moving and extremely slow-moving ground motions, exhibit significant spatial dependence of clustering on active falls.In particular, extremely slow-moving ascending PS are found spatially clustered on active falls at both distances of 18.3-21.3km and 50.8-56.7km.On the other hand, no spatial dependence can be apparently detected between PS and dormant falls in the Arno River basin.In terms of stabilized falls, only very slow-moving descending PS were found evidently clustered on them, with the estimated clustering distances mainly between 77.7 and 80 km.

Creeps
In the landslide inventory, all creeps mapped in the Arno River basin were detected with active states of activities.Therefore, the bivariate K-function analysis was performed on PS and active creeps.The result of the bivariate K-function analysis on PS and creeps is able to reject the null hypothesis that all PSI point targets tend to cluster similarly on active creeps across the Arno River basin (Figure 6, Table 4).In detail, only ascending PSI point targets exhibit spatial dependence on active creeps, with the clustering range of 2.5-16.3km and 12.5-15 km for extremely slow-moving and very slow-moving PS, respectively.No significant spatial dependence can be observed for descending PS on active creeps.

Discussion on Materials
The landslide occurrence data employed in this study were obtained from the landslide inventory produced from multiple sources of datasets.The inventory was mapped at the scale of 1:10,000 throughout the Arno River basin.In particular, temporal aerial photos were employed to facilitate the visual interpretation of landslides.In addition, PS information was added to confirm the states of landslide activities and modify the landslide boundaries.However, uncertainties may be introduced during the preparation of this landslide inventory.These uncertainties may mainly come from the following aspects.
Firstly, with regards to landslides mapped from aerial photos, uncertainties may come from: (1) the quality and scale of the aerial photos and base maps; (2) registration of aerial photos to base maps and the derivative positioning errors of the mapped landslides; (3) the empirical and subjective approach of visual interpretation on the aerial photos for landslide identification and mapping; and (4) the skills and experiences of the involved geomorphologists and geologists during the inventory mapping [40,66,67].
Secondly, in this study the regional topographic maps were used as the base maps on which the landslides were mapped.Therefore, the landslide mapping quality is closely related to the quality of these topographic maps which are expected to provide evident and accurate topographical and morphological features.In this study, the regional topographic maps provide relatively complete and accurate information at 1:10,000 scale especially for recognizing features of roads and buildings, and generally the shape of divides and drainage lines can be clearly identified.However, the changes of land use and the patterns of vegetation may be incomplete due to the fact that the topographic maps used in this study were only updated until 2002.In addition, due to the large scale of mapping (1:10,000), the errors may be further propagated when digitizing the topographic maps and mapping into GIS.
Thirdly, the uncertainties may also been brought by the age and freshness of landslides, and the types of land cover and land use.Also, in general, the landslide inventory does not record all the landslides that occurred, especially for those obscured by erosion, vegetation, urbanization and land use change, thereby making the inventory incomplete and the degree of this incompleteness is often difficult to be judged and quantified [40].
In addition, other uncertainties may be propagated by the PSI technique which at times has limitations in landslide detection and investigation.In particular, due to the 2π-wrapping of interferometric phases accompanied with possible phase unwrapping errors, the measurable velocity between two consecutive images is limited to less than a quarter of the radar wavelength, namely 14 mm in the case of C-band RADARSAT-1 images utilized in this study.For this reason, with regards to the velocity scale of landslides proposed by Curden and Varnes [53], currently the derived PSI point targets may only have usefulness in detecting very slow-moving and extremely slow-moving landslides that are within the detection precision of the PSI technique.This limitation restricts the application of PS only in landslides with very slow dynamics.The major types of landslides in the Arno River basin are slow-moving deep-seated slides (74.8%) and creeps (19.4%).However, rapid-moving shallow landslides and debris flows were also reported in this inventory.As a result, the bivariate K-function analysis on spatial patterns between PS and landslide occurrences may be affected by these uncertainties induced by the detection precision of the PSI technique.Nevertheless, with the new X-band SAR missions (e.g., TerraSAR-X and COSMO-SkyMed) offering higher temporal resolution (3-4 images for each month), the PSI technique can be used to detect higher displacement rates up to 35 cm/year.
Moreover, the PSI technique has limitation in detecting stable benchmarks in areas without stable reflectors and coherence targets, for example, in vegetated and forested areas where natural reflectors are usually unavailable.The Arno River basin is highly urbanized and populated with 2.2 million inhabitants.However, in some landslide-prone areas less affected by human activities, it is sometimes possible to obtain PS with relatively low density.If available, it is suggested to use SAR images acquired from X-band sensors that provide shorter wavelength for the bivariate K-function analysis.It has already been proven that the PS density can be largely increased by processing X-band images such as COSMO-SkyMed [68] and TerraSAR-X data [69].In addition, it is also possible to improve PS density by modifying some algorisms for data processing optimization, such as the SqueeSAR technique that squeezes distributed scatterers for extracting ideal phase vectors [70,71] and multi-master interferograms strategy exploiting sampled coherent scatterer motion at the acquisition time of all available SAR images [69].
Furthermore, it should be noticed that the displacement rates calculated from PS are averages of observation during the time span of the processed SAR images.The process is based on the assumption that the targeted mass movement is linear with constant rate of temporal phase variation [26].Therefore, the temporal variation of ground motion may not be characterized by PSI, which could limit the further understanding of landslide mechanisms.This limitation of PS in describing complex movements of landslides may propagate these uncertainties into the bivariate K-function analysis.

Discussion on Methods and Results
The logical framework behind this study can be structured as: (1) landslides tend to be spatially clustered since they usually occur in places where mass movements have existed in the past [39,40]; (2) landslides can be effectively detected by PSI point targets; and (3) clustering of PS can accordingly provide useful information for landslide detection at certain scale distances.In order to fulfill this concept, the bivariate K-function analysis was utilized in this study, providing a practical tool for analyzing spatial patterns between two groups of point datasets: PSI point targets acquired and processed from RADARSAT-1 images and landslide occurrences mapped in the inventory across the Arno River basin.Spatial heterogeneity was observed among landslides with different states of activities and PS with various orbits and velocities.
In this study, the ascending and descending PSI point targets were analyzed separately in the bivariate K-function analysis.This is because the spatial distribution of PSI point targets may be different when they are processed from various orbits and the bivariate K-function analysis is tenaciously associated with the spatial variation of the analyzed datasets.Alternatively speaking, PSI point targets acquired from one orbit may not be invariably recognized as the same stable benchmarks from another orbit.This can be mainly attributed to two facts.Firstly, the ascending and descending images are processed independently to obtain stable benchmarks and corresponding PSI processing approaches of long-term SAR images for ascending and descending orbits are different.For example, the ascending and descending processing approaches have differences in acquisition dates, master images, reference points, phase unwrapping, coherent points identification and co-registration.Secondly, the SAR processing can often be affected by topographic surfaces.Particularly in mountainous terrains, it is common that radar signals have different behaviors between ascending and descending images.For example, in terms of distortion due to the geometry, it is possible to find that a slope shows strong shadowing effect in an ascending image while the same slope can have significant layover effect in a descending image and vice versa.These effects of layover and foreshortening may influence the recognition of stable benchmarks and accordingly may produce differences in PS location, density and distribution.Hence, it is necessary to separate the ascending and descending PSI point targets for the bivariate K-function analysis.As a result, although the spatial distribution of landslide occurrences remains unchanged, the bivariate K-function rendered different results for the ascending and descending orbits, due to the diverse scatterers acquired and various spatial distributions of those PSI point targets.
In terms of the Ripley's K-function analysis, to achieve the sufficient simulation bandwidth, the high number of Monte Carlo permutations can be suggested [37,48,58].In this study, in order to delimit the envelopes, 99 Monte Carlo simulations of the inhomogeneous random point process were carried out, rendering a confidence level of 98%.Although a more conservative envelope can be achieved for better visualization of predicted curves using the even lower confidence level of 95% (α = 0.05) and 90% (α = 0.1), the corresponding simulation number (39 and 19 times, respectively) of Monte Carlo permutations may be inadequate for a large number of PSI point targets analyzed in the Arno River basin.
It is interesting to note, the results of the bivariate K-function clearly indicated that PSI point targets are not randomly distributed over falls and the clustering behavior was attested over a number of simulation envelopes of Monte Carlo tests.Falls represent a kind of relatively small landslide that removes individual and superficial materials, usually from the cliff surface [72].PSI point targets processed from 4 years of RADARSAT-1 images appear to show spatial patterns of clusterings with active falls throughout the Arno River basin, detected by both ascending and descending orbits with extremely and very slow-moving velocities.This spatial pattern between PS and falls may reflect several factors.In general, the materials detached from the cliff surfaces can hardly be recognized as stable benchmarks by the PSI technique due to the difficulty in identifying targets dropping temporal coherence.In addition to the lack of temporal coherence, the movements of falls are usually not constant and the velocity may exceed the detection precision of PSI.As a result, in the Arno River basin the stable scatterers are mostly identified as extension of outcrops on which the surface displacements may occur within the velocity scales of extremely and very slow-moving movements.Also, the rock outcrops are ideal natural reflectors for the PSI technique that often renders decent scatterers density, thus making them particularly effective in spatial clustering analysis for detection of mass movements at different scale distances.Due to the fact that significant clustering of falls at various scale distances was only detected for active events, the state of activity may have a stronger controlling influence than velocities on the clustering behavior of PS for detecting mass movements of falls.For active falls, extensions of outcrops are prone to be noticed by PS due to the continuous falling processes after the occurrences.In addition, the result shows that the clustering distances for the falls range diversely from small to large scales.It indicates that, in the Arno River basin, the occurrences of falls detected by PS clusterings tend to be connected with controlling factors at diverse scales.In terms of small scales, a group of PSI point targets may be achieved at a localized area, e.g., on a single rock outcrop, thereby bringing an evident clustering of PS.With regard to large scales, the occurrences of falls are likely to be associated with regional hazard mitigation measures, particularly in the Arno River basin, a broad coverage of forests were utilized to protect a large group of inhabitants spread over different communities and provinces, thus making the detection of falls at relatively large scales possible.However, it should be noticed that in this study only six events of stabilized falls were mapped in the landslide inventory and they may not provide sufficient samples for the bivariate K-function analysis.As a result, identification of stabilized falls may be difficult if only using the clustering information of PSI point targets in the Arno River basin.Manual interpretation of PS combined with additional remote sensing data such aerial photos and very high resolution (VHR) satellite images, or field investigation, is fundamentally needed for a more accurate interpretation.
Basically, creeps are prone to be detected by PSI since they represent a kind of slow and gradual shallow mass movement.However, in this study, the results of the bivariate K-function analysis showed that only ascending PSI point targets exhibit significant spatial clusterings with the occurrences of creeps in the Arno River basin.This may result from the topography since in the Arno River basin creeps predominantly occur in the slope aspects oriented in south-west directions and accordingly they are prone to be detected by ascending PS.As a result, spatial clusterings of ascending PS were found evident for creeps in the Arno River basin.This result suggests that the application of PS in landslide studies should not only consider the detection precision but also the acquisition geometry of SAR images and the topographic settings of slopes.In addition, the bivariate K-function analysis suggested that throughout the whole basin those creeps tend to cluster at maximum distances of 16.3 km and 15 km, detected by ascending PS for extremely slow-moving and very slow-moving movements, respectively.Spatial randomness of ascending PS was found statistically evident with creeps above the scale distance of 16.3 km in the Arno River basin.It is considered that some of these clustering effects at relatively small scale distances are possibly associated with the topography of the Arno River basin where high relief mountain slopes are geographically concentrated and the spatial distribution freeze-thaw processes affecting the creeping activities are clustered in space.Moreover, factors controlling the significant spatial heterogeneity of creeps at relatively large scale distances may be associated with climate changes which render ground freezing conditions at the regional scale such as air temperature and precipitation amount [73].
Significant spatial clusterings were detected for descending PS over very slow-moving flows throughout the Arno River basin.Spatial dependence was found evident for PS on active and dormant flows with relative small clustering scales of 1.7-2.7 km and 2.1-4.6 km, respectively.Mackey et al. [74] have summarized the long-term evolution behaviors of flows, seeing that active flows tend to be transient and flow-prone landslide may be featured by historical dormant events.Hence, dormant flows characteristically represent long-term slope evolution in past periods and historical dormant flows may circumscribe a small proportion of active terrain, while active flows may transfer to long-term dormant ones on account of the imbalance between mass entering and exiting the flows [75].Since this phenomenon often exists at a small fraction of terrain, flow occurrences may demonstrate spatial dependence within relatively short scale distances, such as in this study within 4.6 km as detected by descending PS.In addition, it is noticed that descending PSI point targets were found clustered on stabilized flows at relative large scale distances from 73.8 to 78.1 km.However, only 18 events were mapped in the inventory for this type of mass movement, and the scale distances of clustering may be overestimated due to the insufficient data inputs.
Spatial clusterings were also found evident for descending PS on very slow-moving slides throughout the Arno River basin.This is possibly because a large portion of slides in the basin are with slope aspects oriented in east, north-east and south-east directions, which are prone to be detected by descending orbits of the RADARSAT-1 satellite.The significant clustering distances for dormant slides are relatively short within 8.8 km.This could be interpreted as follows.Most of the dormant slides in the Arno River basin occurred at the end of the last glaciations during the early phases of the Holocene as a result of ice retreat [49,76].A large number of these dormant slides are covered by vegetation and accordingly they are very difficult to be recognized as stable benchmarks from temporal SAR images.The vegetated terrains are prevalent in large areas covering those dormant slides and therefore no significant spatial clustering would be detected at a relatively large scale distance.The reason active slides nonetheless were found clustered at relatively large scale distances ranging between 57 and 58.4 km may be ascribed to two factors.Firstly, some dormant slides can be reactivated into active slides due to the intensive rainfall events in the basin and the coverage of these rainfall events may correspond to this scale distance.Secondly, the activities of slides can be altered due to the land use change carrying on a correspondence to this scale distance.In the Arno River basin, one typical example of this land use change is the agricultural activities for wine and olive oil production, which can produce 2.5 billion Euros every year.The number of specialized vineyards and olive groves has greatly increased since the 1960s and they have been developed with agricultural mechanization, extensive slope leveling and alteration of natural drainage network [77].These activities contribute largely to the soil erosion and change the slide activities from dormant to active.No spatial dependence was observed between PS and stabilized slides because mass movements are no longer affected by the original causes and they are below the velocity threshold as defined for extremely slow-moving landslides.
Very slow-moving and extremely slow-moving landslides can be detected through the semi-automatic PS clustering analysis from spatial statistics, which is particularly useful for application in a large area so as to rapidly identify hazardous spots for quick risk mitigation purposes [45,46].For instance, by means of Persistent Scatterers Interferometry Hotspot and Cluster Analysis (PSI-HCA), it is possible to identify those landslide hotspots and further detect an active slide in Trespiano area which was previously unmapped in the inventory (Figure 7).For these feature pattern analyses, the selection of a suitable scale distance is essentially needed.However, how to select this appropriate scale distance remains as a key problem in these studies.For example, in the study of Lu et al. [46], the scale distance was defined as the average value of the shortest paths to a channel and ridge calculated for all pixels.This method may be meaningful for a single slope but is very limited at the regional scale.The bivariate K-function analysis nevertheless is able to investigate spatial patterns over a range of multiple distances from a slope to the regional scale.That is, the bivariate K-function analysis enables a calculation at multiple distances and it reveals how spatial pattern of point distribution varies with scales.Therefore, the bivariate K-function makes possible a selection of the most appropriate scale for spatial pattern analyses such as Moran's I spatial autocorrelation analysis and Getis-Ord G Statistics.In this work, the bivariate K-function analysis on spatial relationships between PS and landslide occurrences is expected to provide a guideline for studies trying to detect and investigate landslide occurrences at regional scale through spatial statistical analysis of PS, from which a priori understanding of the influence of scale distances on the investigated landslides can be achieved.For example, the clustering distances can be firstly examined before performing the landslide hotspot mapping for the study of Bianchini et al. [45].Lu et al. [46] noticed the influence of scales in the PSI-HCA, which utilized the Getis-Ord G Statistics and the kernel density estimation for a rapid detection of extremely slow-moving landslides from the clustering nature of PSI point targets on landslides.However, their study did not differentiate the scale distances on different types of landslides and fails to describe the spatial patterns of PS and landslides.The result of the bivariate Kfunction analysis in this study may provide a reference for determining the scales used in spatial statistical approaches, such as the Getis-Ord G statistics and the kernel density estimation, before the PSI-HCA analysis.In addition, the bivariate K-function analysis may help to understand the usefulness of PSI-HCA on different types of landslides.For example, from the results of bivariate K-function analysis in this study, PSI-HCA tends to be more useful in detecting active falls in the Arno River basin at different scale distances.Therefore, for these spatial clustering analyses of PS targeting on landslide mapping and detection purposes, it is suggested to apply this bivariate K-function in advance so as to test the applicability of these methodologies at the beginning stage.With well-defined scale distances, the clustering behaviors of PS for landslide detection and investigation can be more carefully determined.
Finally, it needs to be noticed that in this study the computation time for the bivariate K-function analysis under the R environment is relatively long using the workstation of a 3.4 GHz processor and a 8G RAM.In this study, the landslide polygons in the inventory were converted to point datasets based on their centroids with one of the motivations to reduce the computation time.If the analysis is based on original landslide polygons, it is computationally difficult to be accomplished with this configuration of workstations.However, this approach may induce uncertainty especially related to the variance in the size of landslides.In the Arno River basin, the bivariate K-function analysis on PS and landslides tends to be less influenced by this conversion considering that the size of landslides is not greatly varied and 99% of landslides in the inventory are within 0.3 km 2 .Nevertheless, in the future, for possible analysis with larger datasets, it is suggested to use computers with higher speed processors and larger memories, or if possible, to employ computer clusters, for the improvement of computation performance.

Conclusion
This study has deployed the Ripley's bivariate K-function to analyze the spatial relationships between four categories of PS having various velocities and orbits (slow-moving ascending, slow-moving descending, extremely slow-moving ascending and extremely slow-moving descending) and four types of landslides (slides, flows, falls and creeps) classified with reference to the classification suggested by Hutchinson [39] across the Arno River basin in Central Italy.The primary purpose is to investigate the spatial patterns between different PSI point targets and landslide occurrences, so as to understand if PSI benchmarks are spatially clustered (spatial attraction) or randomly distributed (spatial independency) with reference to different types of landslides across the basin.In addition, the proposed method is able to assess the degree and scale distance of PS clustering on a variety of landslide occurrences.The results are able to reject the null hypothesis that in the Arno River basin PSI point targets with diverse velocities and orbits appear to cluster similarly with respect to different types of landslides under a 98% confidence interval.In terms of slides, spatial clusterings of very-slowing moving descending PS were detected as statistically significant on active and dormant slides, with clustering scale distances of 57-58.4km and 0-8.8 km, respectively.In regard to occurrences of flows, spatial clusterings were found evident for very-slowing moving descending PS on all types of flows, with detected scale distances of 1.7-2.7 km for active flows, 2.1-4.6 km for dormant flows and 73.8-78.1 km for stabilized flows.As for occurrences of falls in the Arno River basin, all types of PS were found spatially clustered on active falls and randomly distributed on dormant falls, while only descending PSI point targets were observed spatially clustered with corresponding scale distances ranging between 77.7 and 80 km.Finally, it is noticed that significant spatial dependence between ascending PS and creeps exists in the Arno River basin, with clustering scale distances of 2.5-16.3km for extremely slow-moving PS and 12.5-15 km for very slow-moving PS.These results show evident spatial heterogeneities of PS scattered on different types of landslides ranging between various scale distances.Significant influence of PS velocities and their acquisition orbits can be found on landslide processes of triggering, stabilization and reactivation across the Arno River basin.Although the bivariate K-function analysis on spatial relationships between PS and landslides may be influenced by the quality of landslide inventory and SAR images, it is expected to provide a guideline for studies attempting to detect and investigate landslide occurrences at a regional scale through spatial statistical analysis of PS, for which a preliminary understanding of the influence of scale distances on investigated landslides is fundamentally needed in advance.

Figure 1 .
Figure 1.The geographical location and landslide distribution in the Arno River basin.

Figure 2 .
Figure 2. The bivariate K-function investigates spatial patterns of landslides and PSI point targets at multiple distances (h), counting the number of neighboring PS for each landslide within the given range of multiple distances.

Figure 3 .
Figure 3.The results of the bivariate K-function analysis for assessing spatial relationships between PSI point targets and occurrences of slides across the Arno River basin.Dashed lines indicate 98% confidence intervals and solid lines represent the observed values.Observations outside the confidence intervals reveal spatial dependence of clustering (above) or repelling (below) at corresponding distances.Observations inside the confidence intervals indicate spatial randomness.

Figure 4 .
Figure 4.The results of the bivariate K-function analysis for investigating spatial relationships between PS and occurrences of flows in the Arno River basin.Dashed lines delineate 98% confidence intervals and solid lines record the observed values.Observations outside the confidence intervals suggest spatial dependence of clustering (above) or repelling (below) at corresponding distances.Observations inside the confidence intervals denote spatial randomness.

Figure 5 .
Figure 5.The result of the bivariate K-function analysis for analyzing spatial relationships between PS and occurrences of falls within the Arno River basin.Dashed lines express 98% confidence intervals and solid lines signify the observed values.Observations outside the confidence intervals indicate spatial dependence of clustering (above) or repelling (below) at corresponding scale distances.Observations inside the confidence indicate spatial independence.

Figure 6 .
Figure 6.The result of the bivariate K-function analysis for evaluating spatial relationships between PS and occurrences of creeps throughout the Arno River basin.Dashed lines are 98% confidence intervals and solid lines are the observed values.Observations outside the confidence intervals suggest spatial dependence of clustering (above) or repelling (below) at corresponding distances.Observations inside the confidence intervals represent spatial independence.

Figure 7 .
Figure 7.The detection of an active slide in Trespiano through PSI Hotspot and Clustering Analysis (PSI-HCA).(a) The landslide hotspot map covering Pistoia-Prato-Firenze basin and the Trespiano area: blue and red hotspots indicate mass movements moving towards and away from LOS, respectively.(b) The previously mapped slides and a newly detected active slide in Trespiano.

Table 1 .
Datasets used for the landslide inventory mapping within the Arno River basin.

Table 2 .
The types and numbers of mapped landslides within the Arno River basin.

Table 3 .
The utilized images for Persistent Scatterer Interferometry processing

Table 4 .
The results of the bivariate K-function analysis showing spatial clustering or independence between PSI point targets and landslide occurrences.

Moving Very Slow-Moving Extremely Slow-Moving Very Slow-Moving Slides
(+) evidence of significant clustering; (-) no evidence of significant clustering or repelling.Numbers indicate corresponding clustering distances with 98% confidence interval.