Landslides caught on seismic networks and satellite radars

We present a procedure to detect landslide events by analyzing in-sequence data acquired from regional broadband seismic networks and spaceborne radar imagery. The combined used of these techniques is meant to exploit their complementary elements and mitigate their limitations when used singularly. To test the method, we consider a series of six slope failures associated to the Piz Cengalo rock avalanche that recently occurred in the Swiss Alps, a region where we can benefit from high spatial density and quality of seismic data, as well as from the high spatial and temporal resolution of the European Space Agency (ESA) Copernicus Sentinel-1 radar satellites. The operational implementation of the proposed approach, in combination with the future increase in availability of seismic and satellite data, can offer a new and efficient solution to build and/or expand landslide catalogues based on quantitative measurements and, thus, help in hazard assessments and the definition of early warning systems at regional scales.

1 Introduction

Landslides cause fatalities and devastation globally, with remarkable effects especially on low-income and/or developing countries (Froude and Petley, 2018). While the spatial occurrence of landslides is related to intrinsic geomorphological and climatic characteristics (Stead and Wolter, 2015), catastrophic failures arise when slope materials reach a critical damage state (Petley, 2004). In many cases, the ultimate trigger towards failure events is related to anthropic activities, meteorological events, and earthquakes (Bayer et al., 2018; Huang et al., 2017; Lacroix et al., 2019).

Quantitative and accurate data on the timing, location, and size of landslide events are crucial for studying the relationships between local and regional preconditioning factors, to recognize potential causes, and to identify the potential effects of climatic forcing. Moreover, efficient early warning systems at regional scale rely on the availability of accurate and complete landslide catalogues (Gariano and Guzzetti, 2016). Despite recent efforts, the knowledge on spatial and temporal landslide distribution is incomplete. The information about landslide volume, runout, velocity, etc., is usually available only when the events threat life or damage infrastructures, as well as when they are associated with large earthquakes or exceptional meteorological occurrences. These catalogues, however, deliver only a partial picture of the impact of such events on the landscape. In addition, many landslide events are unreported because they occur in remote regions and do not have immediate and/or relevant impacts on human activities. This strongly hinders the completeness of inventories used for hazard assessment and for the calibration of early warning systems at regional scales (Guzzetti et al., 2019).

Figure 1 Overview of the area of investigation. (a) View of Val Bondasca, with an approximate outline of the area affected by the Piz Cengalo (46.29475^∘ N, 9.602056^∘ E) rock avalanche and subsequent debris flows (© Google Earth 2021). (b) Detail of the release area on 25 August 2017. (c) Detail of the deposits on 30 August 2017. Photos: © VBS swisstopo Flugdienst.

In recent years, two methods have emerged in the panorama of landslide event detection, i.e., satellite remote sensing and seismic data analyses. This is mainly due to the increased availability and quality of these datasets at the global scale and to the open data access policies. In particular, Earth observation (EO) data acquired through different satellite missions are more and more exploited by systematic visual interpretation, and supervised and unsupervised automatic classification methodologies, in order to build catalogues of landslide events triggered by large earthquakes and/or extreme meteorological events (Mondini et al., 2019; Tanyaş et al., 2017). These methodologies strongly depend on the availability of the images, which are usually not adequate for systematic early landslide detection. Despite the identification of signatures of landslide events in seismic networks deployed for earthquake monitoring having been occasionally studied in the past (Govi et al., 2002; Weichert et al., 1994), current technical advances and the diffusion of broadband seismic sensors have increased the possibility also to detect and locate landslide events of small–moderate size at regional scales. Automatic or semi-automatic procedures adapted from earthquake location routines have demonstrated fair performances (Chao et al., 2017; Dammeier et al., 2011; Ekström, 2006; Fuchs et al., 2018); however, while uncertainties of several kilometers can be tolerated in case of earthquake epicentral locations, landslides are extremely confined phenomena that affect a single slope (or only small portions of it). A more accurate location of the events can be achieved with local networks specifically designed to identify mass movements (Dietze et al., 2017; Cook and Dietze, 2022). Nevertheless, such procedures are impractical when the areas to explore are very large, and the number of stations is poor, as it is typical at the scale of entire mountain chains.

In this work, we jointly use broadband seismic data and spaceborne radar imagery to show a procedure allowing for a systematic detection and location of landslides, as well as an initial definition of their area of impact, and their magnitude. We present results over the region recently affected by the Piz Cengalo, a steep granitic massive located in the central Alps at the border between Switzerland and Italy (see Fig. 1). The area was repeatedly affected by large ($>\mathrm{1}\times {\mathrm{10}}^{\mathrm{6}}$ m³), rock slope failure processes in the past decades, with the main event on 23 August 2017 being the largest ($>\mathrm{3}\times {\mathrm{10}}^{\mathrm{6}}$ m³) and most catastrophic reported in recent years, causing eight fatalities and damages in the range of EUR 50 million (Andres and Badoux, 2018). A detailed description of the event, its preconditioning factors, potential causes, the dynamics of the rock slope failure, and the subsequent debris flow reaching the village of Bondo, is beyond the scope of this work. Thus, the readers are referred to the recent literature for more information on these specific topics (Mergili et al., 2020; Walter et al., 2019).

2 Materials and methods

We consider Piz Cengalo as an exemplary case to demonstrate the potential and the limits of the combination of seismic and spaceborne radar data to provide quantitative information on landslide occurrence in an alpine scenario. We benefit from the high spatial density of the AlpArray seismic network (Hetényi et al., 2018) and from the unprecedented spatial and temporal resolution of Sentinel-1 Synthetic Aperture Radar (SAR) imagery (Torres et al., 2012).

In the following, we describe the steps to initially define a candidate location region with seismic data and then apply change detection investigations on Sentinel-1 SAR imagery to refine the location and identify the slope failure event. Hereafter, we will use the term “landquake” (LQ) to define landslide events recorded by seismic sensors, as increasingly proposed in literature (Chen et al., 2013). However, this term is not meant to provide additional details on specific landslide characteristics.

Figure 2 Seismic network and data. The left panel shows the AlpArray network of broadband stations (red triangles). The right panels show the selected signals (vertical component HHZ, where H is high broadband, H is high gain, and Z is vertical orientation) recorded by AlpArray stations located at different distances from event LQ2 (see Table 1), which occurred on 23 August 2017 (i.e., the main Piz Cengalo rock avalanche event; yellow star in the map).

2.1 Seismic data processing

We consider a total of six events occurred at Piz Cengalo between 21 August and 10 October 2017. The landquakes are characterized by different volumes and runout and all occurred in the same slope but at different stages of the progressive failure process. LQ1 occurred 2 d before the main failure, with three events on 23 August 2017, (LQ2–LQ4), while LQ5 occurred about a month later and LQ6 about 2 months later). Figure 2 shows the distribution of the AlpArray stations and examples of the signals for the LQ2 detected at different distances from the source. The apparent velocities are of the order of 3 km s⁻¹ and are thus compatible with surface waves generated by surficial mass movements (e.g., Dammeier et al., 2011). The Swiss Seismological Service (SED) routinely recognizes landslide phenomena in the seismic records of stations located in Switzerland and in the vicinity of the national borders. Despite monitoring procedures not being optimized to detect mass movements, these are systematically reported. After an event detection (at least three stations triggered on the SED network), a first-order solution is obtained by (visually) identifying coherent energy at multiple stations, typically due to S waves, and using a regional 3D velocity model to estimate location. In general, locations are more accurate when seismic stations are close to the event and there is good azimuthal distribution of observations. For the Piz Cengalo landquake event associated to the largest failure (LQ2), the closest station recording the event is at ∼25 km, and the location accuracy has uncertainties of the order of ±5 km.

To perform our back analysis on the Piz Cengalo sequence, we arbitrarily define a temporal window of 10 min centered on the date and time provided by SED with the manual procedure described above. We consider the waveforms recorded by all the AlpArray broadband stations available for each event and focused on the HHZ channel (i.e., the vertical velocity component of high broadband sampled at or above 80 Hz, generally 100 or 200 Hz). The choice of the HHZ channel is justified by previous studies showing that such a component usually entails the largest energy in case of landquakes (e.g., Dammeier et al., 2011). We apply a short-term average/long-term average (STA/LTA) detection (see the details and parameters in Table S1 in the Supplement) to find the onset time of the event at each station. Then, we compute the time delay between the first triggered station (i.e., the first station recording an event, which is assumed to be the closest to the event) and all the other stations identifying an event in the same temporal window. The resulting values are interpolated on a regular grid of 0.25×0.25^∘, spatially smoothed with an average filter (3×3 kernel), and then normalized to obtain a new function defined here as the likelihood of landquake location (LLL). The candidate region of interest (ROI) potentially affected by a landquake is defined by considering LLL > 0.95 and constraining the change detection processing on a spatial subset of available Sentinel-1 radar scenes.

2.2 Sentinel-1 SAR data processing

We adopt the change detection approach proposed in Mondini (2017), which is specifically modified here to tackle single events instead of populations of landslides. The analysis is performed to identify potential variations in surface backscattering occurred between the pre- and post-event images over the area with LLL > 0.95. In fact, changes in the radar brightness coefficient (beta nought, β₀) have demonstrated to be a suitable indicator for the detection of landslide events of different sizes and occurring in different geographic scenarios (Mondini et al., 2019). In the maps recording the temporal changes in β₀, landslides appear as clusters of similar values in a bulk of speckles. After data acquisition, pre- and post-event Sentinel-1 radar imagery require the following steps: radiometric and geometric corrections, multi-looking, filtering of the intensity values, co-registration (pair alignment), computation of changes (logarithm of the β₀ ratio), and ellipsoid correction. The result is a β₀ changes map with a cell resolution of about 14×14 m. Furthermore, the β₀ changes map is segmented using a parametric watershed approach (Roerdink and Meijster, 2000) in which the scale level and the moving window kernel size parameters of the intensity algorithm are automatically assigned minimizing a cost function (Mondini, 2017). The segmentation process is aimed at identifying a few large segments (e.g., the largest, potentially delineating changes associated to the landquake) in the candidate area LLL > 0.95 and a number of small segments intercepting the speckle-like effect. Thus, the landquake event can be recognized as an outlier in the segment's distribution of the areas. The boundaries of the outlier segment provide the potential location of the landquake.

3 Results

To distinguish between local earthquakes and landquakes, we applied the method proposed in Manconi et al. (2016), based on the ratio between the local magnitude and the duration magnitude (see Table 1). This method classifies the Piz Cengalo sequence as landquakes. Moreover, the volumes computed by considering the empirical relationship with the duration magnitude (Manconi et al., 2016) are in agreement with the ones measured with lidar (see Walter et al., 2019).

Figure 3 Likelihood of landquake location (LLL) based on the arrival time of seismic signals recorded by AlpArray stations. This basic analysis of the seismic data is used to constrain the approximate location where a landslide event has occurred. (a) LLL over the entire AlpArray network. (b) Magnified view of the areas with high likelihood. The area 0.95 < LLL < 1.0 is used to confine the change detection analysis. The true location of the Piz Cengalo event (white star) is also shown.

Table 1 Summary of the landquakes analyzed in this study and associated to the Piz Cengalo slope failure. ML (local magnitude) is estimated by SED, while average magnitude duration (MD) and volumes are computed, following Manconi et al. (2016), by considering the event duration on all triggered AlpArray stations. Note that all LQ events have ML/MD less or equal to 0.8, i.e., they can be discerned from earthquake events which typically have ML/MD ∼ 1.

Figure 3 shows the results obtained by analyzing the seismic data available for the LQ2 event. This is the largest landquake, and its seismic signature was detected by tens of stations up to ∼400 km distance from the source (also see the Supplement). The computed LLL function is approximately centered on Piz Cengalo massif. The area within LLL > 0.95 is of the order of 35 000 km², i.e., ∼2 % of the entire grid considered in the interpolation. However, this area is still very large for an accurate identification of a slope failure event that affects an area of about 1 km² (Walter et al., 2019). The initial candidate region defined by the LLL function is used to first identify the available Sentinel-1 imagery in terms of time of acquisition and orbit. In this specific case, the suitable Sentinel-1 orbits are the T015, ascending, and T066, descending, respectively. Then, the change detection processing is not applied to the entire image but only to the area with LLL > 0.95, which is 20 % of the acquired SAR scene.

Figure 4 shows the results of the change detection analysis obtained on the ascending T015 imagery (see Table S2). Due to the temporal proximity of the LQ1–LQ4 sequence (which occurred within 2 d; see Table 1), the events cannot be singularly discriminated because the Sentinel-1 constellation (when both Sentinel-1A and 1B are operative) revisit time in Europe was of 6 d in the period of analysis. The LQ2, however, has certainly been the main cause of the surface changes, and for this reason, we refer hereafter mainly to this event. The outlier segment covers an area of ∼0.9 km², about 2 orders of magnitude larger than the average areas of the segment's distribution. Moreover, the segment is elongated and thus has a very low value of roundness (defined as the area of the circle with the same length as the polygon to the polygon area; see results in the scatterplot in the bottom panels of Fig. 4). The footprint and the dimensions of the segment identified are in very good agreement with the area affected by the rock avalanche (Walter et al., 2019). Since the events LQ5 and LQ6 are smaller in magnitude compared to the LQ1–LQ4 sequence, the changes in the SAR image cannot be univocally and automatically defined as for the LQ2 (see Figs. 4 and S2 and S3). In fact, the LQ5 event is likely the fifth in size segment in the region analyzed, with quite a stretched/elongated shape and a roundness value of 0.17 (median roundness equal to 0.47; first quartile equal to 0.37). The first and the third segments are contiguous, rounded in the first and amorphous in the third, indicating a large change, of ∼21.4 km² about 25 km north of Piz Cengalo, in an area where landslides were not reported. The second and the fourth segments, contiguous and rounded as well, record a change of ∼15.3 km² over the Cima di Castello mountain, 6 km east of the Piz Cengalo valley. Their roundness ranges from 1.6 (for the amorphous segment) to 0.45. The LQ6 event does not produce statistically relevant changes in terms of size and in terms of shape.

Figure 4 The top panel shows the results of the change detection analysis. The red polygon shows the area identified as potential landquake location for the main landquake event (i.e., LQ1–LQ4) identified by processing the Sentinel-1 pre- and post-event, while the gray polygon is the area hit by the rock avalanche (compare Walter et al., 2019). The white star and the yellow star show the locations of the largest segments for LQ5 and LQ6, respectively, as identified within the Bondasca valley. Also see the Supplement for more details on the segmentation results. The black dots show the epicentral locations provided by SED (see Table 1). The bottom plot shows scatterplots indicating the distribution of changes in radar backscatter versus normalized areas of the segments. The color map shows the normalized roundness. The segments associated to the landquakes are identified by the black arrows. See the text for more details.

4 Discussion

Seismic data are capable of providing indirect evidence of the time of landslide occurrence, also in inaccessible locations, but independent verification of the location is necessary for event confirmation and classification (Ekström and Stark, 2013). On the other hand, remote sensing data can deliver direct evidence of the areas hit by landslide events, but independent observations are necessary to identify the exact time of occurrence (Guzzetti et al., 2012).

We propose here an approach exploiting seismic and remote sensing (specifically, spaceborne SAR data), which is suitable for the development of automatic pipelines aimed at a systematic identification, location, and first evaluation of landslide events. We have shown, as an exemplary case, the application to a sequence of events that recently occurred in the Swiss Alps. Our results provide several hints on the potential application of this approach in operational scenarios. We have applied a STA/LTA approach for the identification of the event on an arbitrarily constrained temporal window. The STA/LTA method has shown to be suitable for the automatic detection of mass movements in continuous seismic records also for early warning purposes, although specific calibration of the parameters used is necessary and depends on the sensors, the network configuration, and local conditions (Coviello et al., 2019). One of the main arguments against the use of the STA/LTA approach in the detection of mass movement signals lies in the inaccuracy of the determination of the event's onset, which might cause errors in the subsequent location procedures (Fuchs et al., 2018). Since we refine the location using the remote sensing imagery, the STA/LTA approach is sufficient to constrain the candidate region for the change detection task. Inaccuracies up to few seconds of the STA/LTA detection that would cause large inaccuracies in location routines based on seismic data only would cause only negligible changes in the LLL function.

An important problem after the detection of an event is the distinction and/or classification of the signals recorded in continuous seismic waveforms (earthquakes, explosion, mass movements, anthropic sources, etc.). Several authors have proposed empirical-based relationships, signal processing, and/or machine learning strategies, which achieve good performances (Dammeier et al., 2016; Hibert et al., 2014; Moore et al., 2017). The approach used here, proposed by Manconi et al. (2016), shows that, in the Piz Cengalo sequence, all events could have been classified as landquakes. Moreover, the evaluation of the rockslide volumes based on the empirical relationship observed with the duration magnitude provided good agreement with independent volume measurements. The same approach has recently been implemented in an operational regional system in Taiwan and shows encouraging results (Chang et al., 2021).

Despite the candidate location being identified with a basic proximity approach, the source region is already reasonably well constrained for all six LQ events considered (see also Fig. S1 in the Supplement). This result is possible only when a relatively high spatial density of seismic sensors is available, such as the AlpArray network. We tested how the estimation of the area LLL > 0.95 would vary by omitting the nearest seismic stations (see Table 2). For LQ1 and LQ2, there is a clear increase in the LLL > 0.95 area, although the values in percentage with respect to the entire size of the network still show a substantial benefit in the reduction of the candidate area for the subsequent change detection analysis. In case of smaller events, the differences are difficult to evaluate because the number of stations operative at the moment of the event, the signal/noise ratio of the seismic data, and the grid interpolation step have a more relevant effect. The results of this test show that the size of the area for initial guess of landquake location is thus not only related to the distribution of the stations but also to the size (and likely the dynamic) of the event. At this stage, we cannot provide general network requirements that could be exported in other case scenarios. More advanced location routines can be applied, but the homogenization of procedures across large areas like the entire alpine chain is not straightforward. In addition, an increased level of complexity would not certainly correspond to an increase in accuracy for landslide location.

Table 2 Results of the assessment of differences between areas with LLL > 0.95 when reducing the number of stations used for the interpolation. In case 1, we removed all stations triggered between 0 and 5 s, i.e., up to ∼50 km from the event. In case 2, we removed stations triggered between 0 and 10 s, i.e., up to ∼65 km. Missing values are because of lack of stations for reliable interpolation. Values in percentage are computed with respect to the entire interpolation grid, i.e., 1.71×10⁶ km².

Table 3 Summary of the pros and cons related to the herein proposed strategy for landquake identification, location, and classification.

As far as the change detection analysis of the Sentinel-1 SAR data is concerned, the location of landquakes in LQ2 (i.e. in this case the LQ1–LQ4 sequence) is straightforward. The event was large and caused a vast drop in the backscattering coefficient in the post-event image, which occurred by spatially oversizing the surrounding random changes always present in SAR images (speckle-like effect). Furthermore, other environmental changes in the area are not relevant and, in this specific case, are mostly in the direction of an increase in the backscattering coefficient. The results of the segmentation are unambiguous in all the images, whatever the acquisition mode and the polarization are, even if the final segments can be slightly different. Additionally, post-processing, like smoothing or gap-filling filtering, can also partially change the final shape of the segment and the identified area. On the contrary, the identification of the LQ5 and LQ6 events shows more complexity, and it is not straightforward. Their signals emerge only in the ascending imagery with VH (vertical–horizontal) polarization, which is a possible indication of a weak change in roughness along the slope (Sung and Holzer, 1976). According to seismic data, their sizes are much smaller compared to LQ2, and then the corresponding changes of the backscattering coefficient are expected to be less distinguishable in the bulk of random speckles (see Fig. 4). In fact, when the signs left in the β₀ changes map have a size comparable to other environmental changes or the speckle-like segments, landslides cannot be univocally recognized. Regarding LQ5, only a supervised post-processing, including customized filtering to facilitate the segmentation, and manual cluster shape analysis (Mondini et al., 2019) over the valley allowed the highlighting of a potential segment of interest. The segment is the fifth in size among millions, with a stretched shape compatible with the slope process under study. On the contrary, the first four clusters are rounded, or amorphous, and more adequate for representing other types of processes that occurred in the area. Multi-band, multi-polarization data, and a shorter revisit time of the satellite would have probably helped in reducing the environmental noise and then in surfacing the segment. For LQ6, a small but clear signal along the slope is present in the catchment, but it is not large when considering the entire distribution of segments. Other geometrical parameters such as, for example, the elongation or the roundness of the segment's area do not help (see Fig. 4). We can then assert that LQ6 is below the limit of the spatial resolution of the used images. A potential adaption for the operational implementation of our approach could be running the change detection task on progressively increasing LLL thresholds (0.95, 0.975, etc.). This could provide additional hints on possible hot spots, which can be verified with subsequent SAR acquisitions and/or supplementary remote sensing imagery (spaceborne or airborne).

In Table 3, we summarize the strong points and the shortcomings of the herein proposed methodology, including the working hypotheses to overcome the current limitations. One of the main advantages of the pipeline is that the data pre-processing can be fully automated. Moreover, some of the limitations depend mainly on technical constraints that might be overcome in the near future. There are a number of problems, however, that depend on the intrinsic limitations of the data considered, and it is difficult, at the current stage, to have a clear definition of the best-suited strategies for improving the performance. We will perform future evaluations on the continuous processing of seismic data to make a substantial assessment of the potential implementation of our procedure, to find the best compromise for specific parameters, and to minimize false/missed alarms.

5 Conclusions

The key take-home message of this study is to show how the systematic combination of seismic and remote sensing data can be useful for identification and mapping of landslide events. The use of SAR satellites shows the advantages of all weather, day and night, and systematic acquisitions at the global scale. When available, optical imagery and/or SAR imagery acquired with different bands, fully polarimetric, or with a higher spatial resolution can eventually contribute to increase the quality and the quantity of the information. We believe that the combination of seismic and spaceborne data is a viable approach for a future operational monitoring system at the scale of the Alps, and for this reason, this work can be the starting point for raising awareness in the community and to foster the cooperation and funding necessary for such an endeavor. We conclude by remarking that our approach is intended to be used for systematically populated landslide catalogues relying on quantitative and accurate information on timing, magnitude, and frequency and also in remote areas. Improved catalogue completeness is very important for the calibration of regional early warning systems based on rainfall thresholds and on regional hazard assessments (Guzzetti et al., 2019). An increase in the availability of remote sensing imagery with daily or sub-daily revisit times could lead to an employment in early detection of landslide events and possibly in disaster response scenarios, but these potential applications have to be evaluated in future studies.