SAR Data for Detecting Landslide Phenomena: The November 26, 2022 Landslide of the Ischia Island (Southern Italy)

In this article, several change detection techniques based on satellite SAR data acquired by Cosmo-SkyMed Second Generation missions have been evaluated aiming at detecting the landslide-mudflow phenomenon triggered by the strong flooding event that hit Ischia Island in November 2022. It severely impacted on Casamicciola Terme area causing damages and collapses of many buildings and unfortunately also casualties. Experimental results show how both single-polarimetric (SP) and dual-polarimetric (DP) techniques are able to detect the main landslide occurring along the northern flank of Mt. Epomeo and in some cases also connected phenomena such as mud accumulation. Additional analyses have been performed to quantitatively evaluate the performance of all of them showing, as in this case, DP techniques outperform the SP ones. The outcomes are then discussed taking into account both the features of each technique and the investigated scenario. Detecting and mapping this kind of phenomenon is important for the evaluation of the affected area, especially for complex scenarios such as Ischia island, and can be very useful to support both the stakeholders for first aid and the civil protection for the post-crisis management.

SAR Data for Detecting Landslide Phenomena: The November 26, 2022 Landslide of the Ischia Island (Southern Italy) I. INTRODUCTION O N NOVEMBER 26th, strong flooding hit Ischia island, in the Gulf of Naples (southern Italy), causing several problems to the infrastructures and triggering a landslide characterized by huge mudflow on Mt. Epomeo, in the hamlet of Casamicciola Terme (see Fig. 1). Unfortunately, the mudflow violently hit the urban area of Casamicciola Terme, completely destroying some buildings, severely damaging and submerging more than thirty private residences, and dragging several cars into the sea. The search for the missing by the Italian Civil Manuscript [12]. The brown polygon indicates the mapping of the damaged areas (grading product) due to the flooding, retrieved by the Copernicus Emergency Management Service (https://emergency.copernicus.eu/). The yellow star represents the epicenter of the earthquake occurred on August 21st, 2017 from INGV source (Osservatorio Nazionale Terremoti). The black lines are the faults from [13].
Protection, the Fire Brigade, and volunteers lasted several days, making it even more difficult due to adverse weather conditions and unfortunately, in the end, there were 12 casualties. The Italian government promptly declared a state of emergency and allocated 4 million euros to provide immediate relief to the local population affected by the landslide-mudflow with more than 200 residents left homeless. It is not the first time that Ischia has been at the center of media attention, being it an island of volcanic origin characterized by ground deformation and seismic activity which are constantly monitored by seismic and geodetic ground-based networks of the Istituto Nazionale di Geofisica e Vulcanologia (INGV) [1], [2], [3] and through remote sensing data [4]. In August 2017, many damages and collapses were registered in Casamicciola following a shallow M w 3.9 earthquake, which also triggered another huge landslide on the same flank of Mt. Epomeo [5], [6], [7]. Now, the landslidemudflow interested a more localized sector of Mt. Epomeo, indeed it was channeled into a narrow drain, likely formed due to the morphology and previous water erosion phenomena, before reaching the urban area of Casamicciola via a road called Via Celario. Mapping the path of the landslide-mudflow is not straightforward but of primary importance to better understand the dynamic and the causes of the phenomenon and to prevent future natural disasters. This is generally performed by visual inspection of optical satellite imagery although they are acquired only with sunlight and can be largely affected by cloud coverage [8], [9]. In this context, satellite data acquired from synthetic aperture radar (SAR) microwave sensors can provide a useful instrument offering a synoptic view of the observed areas sunlight independently, in any atmospheric conditions, and with a temporal sampling of a few days. Several works on the use of SAR data for landslides detection were already performed; however, they focused on the detection of huge single events or the classification of large areas interested in landslide phenomena [10], [11]. Instead, the Ischia case study is quite challenging since the natural channel where the landslidemudflow triggers is narrow (a few meters) and hard to be detected considering the resolution of the SAR sensors. Moreover, the investigated scenario is very complex being characterized by rough topography, vegetation, mud, and urban areas. Then, in this work, several SAR-based techniques applied to X-band Cosmo-SkyMed Second Generation (CSG) images have been investigated to explore the capability of the satellite radar signal to map and detect such kinds of phenomena. Depending on the revisit time of the acquisitions, such kind of information can be exploited for rapid mapping and hazard assessment purposes thus supporting the stakeholders in the rescue of missing and in the management of postcrisis scenarios.

II. GEOLOGICAL SETTINGS
Ischia is an active volcanic island located in the Gulf of Naples (Southern Italy) belonging to the Phlegraean Volcanic District, which also includes the Campi Flegrei caldera, Procida, and Vivara islands. The Ischia island represents the subaerial portion of a volcanic complex characterized by a long eruptive history since at least 150 thousand years Before Present (ka B.P.) [14]. The main activity has been the Green Tuff eruption (55 ka B.P.) with the consequent collapse of the caldera in the central sector of the island. After this event, a large resurgence process resulted in the emplacement of the Mt. Epomeo block starting from 30 ka B.P. [15], with the maximum uplift in the order of about 900 m. The resurgent block is bordered by a normal faults' system producing a morphologically well-defined Holocene graben structure at the base of the northern flank of Mt. Epomeo [16]. Despite this large resurgence interesting the central part of the island, the south-central part is characterized by general subsidence pointed out by geodetic data [17], [18]. The last period of the volcanic activity started 10 ka B.P. with 46 different eruptions mainly between 2900 years B.P. and 1302 A.D. (the Arso eruption). As a consequence, Mt. Epomeo's flanks are characterized by slope instability at different scales from shallow mass movements to large rock and debris avalanches. These phenomena typically occurred during eruptions, earthquakes (e.g., during the last seismic sequence at the end of August 2017), or extreme meteorological events as it happened on November 26th, 2022 [4], [5], [6], [7], [19], [20].

III. SAR DATA
Because of the small size of the phenomenon to be investigated, the high-resolution X-band SAR data are the most suitable to be applied in the present study. Therefore, the acquisitions from CSG mission were exploited since they offer a pixel ground resolution of about 2 × 2 m, which is one of the best in the framework of free accessible non-commercial products. The exploited dataset consists of five dual-polarimetric (DP) (HH+HV) Single Look Complex (SLC) images acquired in strip map mode between October 30th, 2022 and January 1st, 2023 along the descending track and with a revisit time of 16 days. The geometry of acquisition of all the images is characterized by an incidence angle of about 27°and an azimuth angle of about 11°. To reduce the speckle noise characterizing the SLC SAR images, despeckling filtering [21], [22] using a window of 9 × 9 pixels along both range and azimuth was preliminary applied. In Fig. 2, the chart of both perpendicular and temporal baselines between all the pairs is shown. According to the scheduling of the CSG mission, there is an acquisition every 16 days, whereas the orbital tube of the mission is not so narrow with values of perpendicular baseline up to 1000 m. For geocoding the output products, the Digital Elevation Model (DEM) by Nappi et al. [12] was used. It has been obtained by digitizing contour maps and spot elevation in vector format (scale 1:2000) provided by the Consorzio Intercomunale Servizi Ischia. In particular, the triangulation irregular network interpolation technique through ArcGIS ESRI Software was used. Vector data of the coastal line were derived from both the digital color orthophoto (pixel 1 m) of Compagnia Generale Riprese Aeree, airborne 2000, and an IKONOS satellite image (RGB, pixel 1 m) acquired on August 9th, 2000. This process produced a 2 × 2 m DEM with an elevation range spanning 0-787.4 m. The DEM is referenced to UTM zone 33 with the WGS-84 ellipsoid.

IV. DETECTION METHODS
In this section, the theoretical rationale that underpins the proposed multipolarization methodologies aiming at imaging any changes due to the landslide is described.

A. Coherence Variation
The coherence maps between two SAR images indicate the correlation degree of the images. They can be estimated, pixel by pixel, according to the following equation: where x 1 (k) and x 2 (k)are the complex values of the pixel k for the first and the second images, x * is the complex conjugate of x, and < > is the averaging module. The coherence is a parameter sensitive to any changes occurring in the investigated area along the time interval between the acquisitions, impacting both the phase and the amplitude of the SAR signal. It is a normalized metric with values spanning between 0 and 1. High values of the coherence represent quite unchanged scenarios according to the incidence SAR signals, indeed they are typically found in urban areas. Conversely, lower values are typically related to scenarios that change very quickly in time producing temporal decorrelation effects such as sea or vegetated areas. Then, modifications of the observed scenarios due to natural or anthropogenic phenomena such as the landslides-mudflow that hit Ischia island should theoretically modify the electromagnetic response of the soil producing some changes in the coherence values. Depending on the features of the investigated scenario, the detection methods based on the coherence are usually performed by using three SAR images: 1) two images pre-event and one image postevent to evaluate the coherence loss during an event [23], [24] or 2) one image pre-event and two images postevent to evaluate any coherence increase after the event [10]. Here, the second approach was conceptually followed by using four SAR images, a pair before and a pair after the event. Indeed, the co-event pair was discarded due to the large perpendicular baseline of about 600 m, which could affect the results. Then, the coherence pre-and post-event were estimated and any variations of such parameters can be highlighted by considering the difference of the two:

B. Differences in the Radar Backscattering
An electromagnetic wave that travels through time and space can interact with a specific target. During this interaction, some of the energy carried by the incident wave (E I ) is absorbed by the target while the remaining energy is emitted as a scattered electromagnetic wave (E S ) [25].
One way to quantify this interaction is using the normalized radar cross section (NRCS), which measures the amount of electromagnetic energy scattered back to the radar by the target. The NRCS depends on several factors including the wave frequency and polarization. Hence, to take into account polarization effects, if we denote by p the polarization of the incident field (E I ) and by q the polarization of the scattered field (E S ), the polarization-dependent NRCS can be defined as follows [25]: where {p, q} ࢠ {H, V} when a linear horizontal (H)-vertical (V) orthogonal basis is adopted, d represents the distance between the SAR and the target, A 0 is the illuminated area, and | | is the modulus operator. NRCS, also called the backscattering coefficient, is highly sensitive to the type of target included in the observed scene [25]. This means that different targets or targets modified over time due to natural or anthropogenic effects give different backscattering values. Following this rationale, the change detection approach aiming at evaluating any changes over the observed scene due to the landslide-mudflow affecting Ischia island is based on the differences in the polarized backscattering intensity where the subscripts pre and post are related to imagery collected before and after the event, respectively, and I stands for the NRCS evaluated for the copolarized (HH) and cross-polarized (HV) channel Note that, (4) is normalized with the sum of the features evaluated before and after the event.

C. Difference of Covariance Matrices
The observation of changes induced by natural or anthropogenic phenomena including landslides and flooding can be also performed using a DP scattering-based change detection approach proposed in [26] and [27]. In such case, a covariance change detection matrix based on the difference between two covariance matrices collected before and after the event is constructed [26] where C 1 and C 2 are the covariance matrices that refer to DP SAR scenes collected before and after the event, respectively. To find the scattering mechanisms ω − that maximize the changes from the pre-to the post-event imagery, a Lagrange optimization procedure is undertaken that results in the following eigenvalue problem [28]: where λ j (with j = 1, 2) are the eigenvalues that maximize/minimize the difference between two covariance matrices. The largest eigenvalue λ 1 of the matrix C CD is mainly positive, which means that it represents the power of scattering contributions that increase or appear in the scene observed in the post-event, whereas the smallest eigenvalue λ 2 is mostly negative, which means that it represents the power of scattering mechanisms that are reduced or disappear in the scene observed in the post-event [26]. A complex event, such as a landslide, together with the geometry of acquisition of the SAR, needs complementary information given by both λ 1 and λ 2 to better detect the fracture shape due to the landslide. Hence, to take into account every possible change, the two eigenvalues are properly combined by arithmetic mean given by [26]

D. Reflection Symmetry
The second DP change detection approach exploits the property of reflection symmetry [29], and it is based on the amplitude of the interchannel correlation given by where S xx and S xy are the co-polarized and cross-polarized scattering elements, respectively.
This property is verified when dealing with naturally distributed scenes, such as vegetated areas. In fact, over this kind of scenario, a very low correlation between co-polarized and cross-polarized channels is expected, i.e., r about 0 [30]. On the other hand, when a phenomenon like a landslide is observed, it is expected that it breaks the natural scenario and, as a consequence, reflection symmetry no longer applies. This implies that non negligible r values are expected [30]. For the purpose of this study, r is expected to be close to zero when dealing with a natural scenario, whereas significant departures from this reference condition are expected over areas where landslides occurred. Hence, to detect the area affected by the landslide-mudflow, the following change detection feature is considered [31]: where r pre and r post are the reflection symmetry maps evaluated before and after the event. This suggests large (low) Δr values over areas affected (not affected) by natural or anthropogenic phenomena.

A. Multilooked Intensity (MLI) Images
A first prompt analysis consists of a simple visual comparison between the multilooked intensity (MLI) of the images acquired just before and after the event. Then, in this case, the CSG SAR images acquired on November 15th and December 1st, 2022 were considered. Indeed, any changes in the investigated area affect the backscattering coefficient thus producing different intensity signals. Such visual comparison is shown at the scale of Ischia island (see Fig. 3) and in a region of interest (ROI) focusing on the area affected by the event (see Fig. 4). Different intensity values are observed along the flank of Mt. Epomeo interested by the landslide-mudflow. In particular, a section of the path followed by the phenomenon, also evidenced by Copernicus Emergency Management Service  (https://emergency.copernicus.eu/), is quite clear and it is highlighted by the red ellipse in Fig. 4.

B. Coherence Variation
Any changes in the area affected by the event, due to the landslide, the mud covering the urban area of Casamicciola, and to some buildings collapse, should produce a coherence variation. Then, to perform such analysis, two pairs of CSG SAR images, acquired before and after the flooding event of the 26th, were exploited considering the difference of the coherence post-and pre-event.
In particular, the two pairs considered for the coherence analysis consist of the images acquired on October 30th and November 15th, 2022 for the pre-event and on December 17th, 2022 and on January 2nd, 2023 for the post-event. Indeed, the same temporal baseline of 16 days between images should guarantee the same effects of temporal decorrelation for the two pairs. Instead, the pair co-event, i.e., the one formed by the images acquired on November 15th and December 1st, and the pair just following the event, i.e., the one formed by the images acquired on December 1st and 17th, were discarded due to the huge perpendicular baseline value of about 600 and 1000 m, which could introduce some geometric decorrelation effects thus affecting the results. The result of such analysis is shown in Fig. 5. As can be noticed, there is an increase in the post-event coherence which seems to define quite well the path of the landslide-mudflow along the flank of Mt. Epomeo. Moreover, some localized patterns of increased coherence are also observed in the surroundings, especially considering the urban area of Casamicciola located in the northern part of Mt. Epomeo.

C. Normalized Radar Cross Section (NRCS)
The backscattering analysis has been performed with the CSG images acquired on November 15th and December 1st using both the co-polarized (HH) and cross-polarized (HV) channels. The metric ΔI is depicted in Fig. 6(a) and (b) for the HH and HV channels, respectively. By visually inspecting, it can be noted that both channels detect the changes associated with the landslide-mudflow. Although the results are pretty similar to each other, little differences between co-polarized and crosspolarized channels are visible. The pattern is well emphasized with respect to the surroundings and it is represented as a red and a blue stripe for positive and negative values, respectively. The red stripe is related to an increase in the backscattering signal, whereas the blue stripe represents backscattering values lower in the post-event image than in the pre-event one. Such behavior is also consistent with both the reference black polygon and the MLI images shown in Figs. 3 and 4 where, from west to east, two stripes of increased and decreased intensity values are detected in the image acquired on December 1st, 2022.

D. Differences of Covariance Matrices
The first DP analysis has been performed with the CSG images acquired on November 15th and December 1st. The λ image False-color λ imagery obtained from the CSG dataset collected on November 15th and December 1st. The black polygon indicates the mapping of impact and spatial extent of the landslide-mudflow (delineation product) retrieved by the Copernicus Emergency Management Service (https://emergency. copernicus.eu/). (9) is depicted in false color in Fig. 7, where the area affected by the landslide-mudflow is enclosed in the black polygon. By visually inspecting, it is possible to observe how the pattern is well emphasized with respect to the surroundings and it quite well matches the reference black polygon. The detected pattern is represented as a red stripe and a blue stripe for positive and negative values, respectively. The red (blue) stripe is related to an increase (decrease) in the power of the scattering mechanisms over the post-event observation. Even in this case, such behavior is also consistent with the MLI images shown in Figs. 3 and 4 where, from west to east, two stripes of increased and decreased intensity values are detected in the image acquired on December 1st, 2022.

E. Reflection Symmetry
The second DP analysis has been performed with the CSG images acquired on November 15th and December 1st. The change detection method Δr (11) is shown in Fig. 8 where the area affected by the landslide-mudflow is enclosed in the black polygon. By visually inspecting, it can be noted that the pattern is well distinguished from the surroundings, and it quite well matches the reference black polygon. The detected pattern is represented as red (negative value) and blue (positive value) stripes, which correspond to an increase and decrease of the symmetry over the post-event observation, respectively. This result can be explained by the fact that the landslide-mudflow occurred over a natural scenario where the reflection symmetry property is satisfied (values of Δr ≈ 0). Hence, under this condition, such a phenomenon breaks the symmetry of the scene resulting in Δr values far from 0.
Once again, according to the previous analysis, such behavior is also consistent with the MLI images shown in Figs. 3 and 4 where, from west to east, two stripes of increased and decreased intensity values are detected in the image acquired on December 1st, 2022.

F. Quantitative Analysis
To evaluate the performance of multi-polarimetric features in landslide-mudflow detection, a target-to-clutter ratio (TCR) analysis is addressed, which is defined as follows [32]: where X ∈ {ΔI HH , ΔI HV , Δγ, λ, Δr} is the mean value of the feature evaluated within the ROI excerpted over the landslide (subscript t) and the surrounding environment (subscript c).
Results, depicted in Fig. 9, show that the DP features, i.e., λ and Δr, outperform the single-polarimetric (SP) ones in terms of landslide/surrounding environment separability, with TCR values up to 4.5 dB. In particular, the reflection symmetry technique based on the change detection feature Δr shows the best performance with a separability around 1 dB more than the best SP ones represented by ΔI HH . The worst performance is shown by the coherence-based technique, Δγ, with a TCR value around 2 dB less than Δr. It can be also noted that the copolarized feature (ΔI HH ) slightly outperforms the cross-polarized one (ΔI HV ).
To deeper analyze the performance of the detection methods, their behavior is evaluated along the two transects shown in Fig. 10 crossing approximately from north to south (NS) and from east to west (EW) the landslide/surrounding environment boundary.
Both the NS and the EW transects show how ΔI HH , ΔI HV , λ, and Δr are characterized by two peaks with respect to the background corresponding to the two stripes of increased and decreased values related to the landslide. In addition, these techniques show the same trend along both transects with similar values of bandwidth of the two peaks representing the landslidemudflow. However, in terms of amplitude, the DP techniques, Δr and λ, outperform up to about 50% of the SP ones with the best performances shown once again by Δr. On the other hand, Δγ shows only one peak corresponding to the stripe of increased coherence after the event. The amplitude is comparable with the one of Δr but it is characterized by a greater bandwidth especially considering the EW transect where it is about twice that of the other techniques. This means that a larger number of pixels are sensitive to variations of coherence instead of the other parameters.

VI. DISCUSSION
SAR-based techniques adopted in the present work show their capabilities in detecting the landslide-mudflow causing several damages on Ischia island. All the retrieved outcomes are consistent with the rapid mapping analysis performed by the Copernicus Emergency group.
The pattern retrieved by Δi HH , Δi HV , λ, and Δr analysis is quite similar. They all are able to detect the upper part of the phenomenon along the flank of Mt. Epomeo showing, from East to West, two stripes of increased and reduced values of intensity or symmetry according to the considered techniques. On the other hand, the coherence analysis shows an increase of such parameters in the post-event image pairs almost corresponding to the stripe of increased values of the other techniques. Moreover, it also allows the detection of some effects of the phenomenon mainly due to the presence of mud, in the lower part of Mt. Epomeo, closer to the Casamicciola urban area. Such behaviors can be explained by taking into account both the geometry of view of the CSG SAR sensor according to the topography of the area and the signal scattering properties depending on the features of the investigated scenario.
The presence of large vegetated areas before the event returns in a mechanism of backscattering known as volume scattering [33] characterized by quite higher values than the scattering due to a rough surface such as the soil composed of mud, water, and debris after the event. Moreover, the landslide-mudflow activated along a channel where the topography shows an incision of several meters previously formed likely due to water erosive effects. Then, after the event which destroyed the vegetation, such morphology should produce geometric distortions typical of SAR images known as shadowing, layover, and foreshortening [34], because of the incidence angle of about 27°and the geometry of view observing from the east to west along the descending track of the CSG SAR sensor. In particular, by observing Fig. 11, it is clear how the eastern section of the channel undergoes a shadowing effect with typically low values of backscattered intensity and symmetry. Instead, the eastern section is mainly affected by foreshortening and layover effects which produce a bright stripe due to the sum of several contributions in the same backscattered pulse with consequent high values of intensity and symmetry in the post-event image. Then, the combination of the lack of vegetation and the SAR geometric distortions produces the pattern detected by Δi HH , Δi HV , λ, and Δr techniques with two stripes of increased and decreased values.
In a similar way, the increased values of coherence along the western side of the channel are strictly connected to the vegetation. Large vegetated areas such as the one before the event are typically subjected to temporal decorrelation effects which strongly affect the SAR signal causing coherence loss. Conversely, in the post-event image, the coherence increases since the vegetation was completely destroyed in the channel by the landslide-mudflow and the soil does not change its electromagnetic response to the SAR signal. The same effect is observed in the lower part of the phenomenon however the scenario is quite different since that area was characterized by some agricultural fields and residential houses with gardens before the event which still cause coherence loss and a large presence of mud post-event which was able to agglomerate due to the lower slope and produce a coherence increase.
Finally, the coherence increase shown in the surrounding is likely due to the different perpendicular baseline characterizing the two pairs, which is about 250 m and 57 m for the pre-and post-event pair, respectively.

VII. CONCLUSION
In this work, several data processing techniques were applied to a dataset of SAR images acquired by CSG missions to detect the landslide-mudflow event that hit Ischia island on November 26th. The aim has been to fully explore the capabilities of both the amplitude and the phase of the SAR signal to map any changes on the ground induced by the event. Both SP and DP configurations have been exploited returning five detection methods with several features. Experimental results show how all the techniques are able to detect the main landslide event with the best performance highlighted by the one based on the reflection symmetry property. On the other hand, accumulations of mud in areas characterized by a low slope are better constrained by the one based on coherence variation. This study demonstrates how, depending on the geometry of the view of satellites and the features of the investigated scenario, a detection method can work better than the others. Then, in such complex cases, it is advisable to use all the techniques to have a synoptic view of the ongoing phenomenon as much reliably as possible. This study also demonstrates how useful can be satellite data for emergency support in case of natural disasters. The time window ideally indicated to provide useful information for the management of postcrisis scenarios, such as the mapping of Authorized licensed use limited to the terms of the applicable license agreement with IEEE. Restrictions apply. the areas affected by landslide phenomena, is within 2 weeks from the event [35]. However, in some cases, satellite revisit time in the order of 15-16 days could not be appropriate, especially considering X-band data and complex scenarios as mountainous areas strongly affected by tropospheric artifacts and decorrelation effects due to vegetation or snow.
In this sense, it is important to remark the importance of the availability of increased satellite data both for the present and future space missions at least for those areas potentially affected by natural hazards such as Ischia island.