Next Article in Journal
Validation of Sentinel-3 OLCI Integrated Water Vapor Products Using Regional GNSS Measurements in Crete, Greece
Next Article in Special Issue
Mapping Post-Earthquake Landslide Susceptibility: A U-Net Like Approach
Previous Article in Journal
Management Zone Delineation for Site-Specific Fertilization in Rice Crop Using Multi-Temporal RapidEye Imagery
Previous Article in Special Issue
Temporal-Spatial Frequency Rockfall Data from Open-Pit Highwalls Using a Low-Cost Monitoring System
 
 
Search for Articles:
Title / Keyword
Author / Affiliation / Email
Journal
Article Type
 
 
Advanced Search
 
Section
Special Issue
Volume
Issue
Number
Page
 
Journals
Remote Sensing
Volume 12
Issue 16
10.3390/rs12162605
remotesensing-logo
Submit to this Journal Review for this Journal Propose a Special Issue
► Article Menu

Article Menu

share Share announcement Help format_quote Cite question_answer Discuss in SciProfiles thumb_up
...
Endorse textsms
...
Comment
first_page
settings
Order Article Reprints
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Integration of Digital Image Correlation of Sentinel-2 Data and Continuous GNSS for Long-Term Slope Movements Monitoring in Moderately Rapid Landslides

by
Marco Mulas
1,*,
Giuseppe Ciccarese
1,
Giovanni Truffelli
2 and
Alessandro Corsini
1
1
Department of Chemical and Geological Sciences, University of Modena and Reggio Emilia, via Giuseppe Campi 103, 41125 Modena, Italy
2
Emilia-Romagna Region—Regional Agency for Civil Protection and Territorial Security, Strada Giuseppe Garibaldi 75, 43121 Parma, Italy
*
Author to whom correspondence should be addressed.
Remote Sens. 2020, 12(16), 2605; https://doi.org/10.3390/rs12162605
Submission received: 6 July 2020 / Revised: 7 August 2020 / Accepted: 10 August 2020 / Published: 12 August 2020
(This article belongs to the Special Issue Slope Stability Monitoring and Investigation Using Remote Sensing Techniques)
Browse Figures
Versions Notes

Abstract

:
This work explores the advantages and drawbacks of the application of Digital Image Correlation (DIC) to Sentinel-2 Multi Spectral Instrument (MSI) data in conjunction with continuous Global Navigation Satellite System (GNSS) monitoring. The goal is to retrieve a spatially distributed and long-term time-series of slope movements in large-scale moderately rapid landslides. The short revisit time of Sentinel-2 satellites (5 days since March 2017 and 10 days before) increases the availability of cloud and snow free satellite acquisitions of the area of interest, which is a prerequisite for the extrapolation of slope movement time-series using DIC techniques. Despite the Sentinel-2 limited spatial resolution, the derived long time-series can be integrated with—and validated by—continuous GNSS monitoring data. This allows to effectively monitor landslide movements that are too fast for the application of interferometric approaches. In this study, we used the Normalized Cross Correlation (NCC) digital image correlation technique by 51 Sentinel-2 MSI scenes (band 4 with 10 m spatial resolution), acquired between 19 February 2016 and 16 July 2019, to derive the slope movement time-series of the Ca’ Lita earthslide-earthflow in the northern Apennines (Italy). During the period considered, the landslide experienced two to three months-long phases of moderately rapid velocity (around 10 m/month) and, in between, prolonged periods of slow movements (approx. 10 cm/month). NCC results have been integrated with, and are compared to, time series from three continuous GNSS devices located in different geomorphic zones of the landslide. On this basis, the errors and limitations associated to NCC time series are analysed and discussed together with their advantages and potentialities for assessing the spatial distribution and monitoring slope movements during moderately rapid reactivation events.

1. Introduction

A wide range of techniques is nowadays available for monitoring slope movements: (i) in situ instrumentation such as inclinometers, wire extensometers, crackmeters, seismic sensors and Global Navigation Satellite Systems (GNSS) [1,2,3,4,5]; (ii) proximal sensing systems such as Terrestrial Laser Scanner (TLS), Ground-Based Synthetic Aperture Radar (GB-SAR) and Automated Total Station (ATS), which allow us to obtain spatialized data using a single acquiring sensor [6,7,8,9,10]; (iii) unmanned aerial vehicles (UAV) that are able to deliver spatially distributed information at a relatively low cost [11,12,13,14,15]; (iv) remote sensing techniques based on satellite data that, by using specific processing methods (radargrammetry, multi-temporal interferometry and digital image correlation), can deliver ground displacements time-series over large areas [16,17,18,19]. Some of these techniques can be effectively adopted to monitor moderate to rapid moving landslides. Specifically, proximal and remote sensing techniques can be optimal tools to retrieve long-term displacement time-series in inaccessible sites, as they allow us to avoid installation and maintenance of in situ systems [20,21].
Among remote sensing techniques, Digital Image Correlation (DIC) algorithms are quite versatile, since they can be used with various kinds of multitemporal georeferenced datasets. This class of algorithms (also known as Offset Tracking—OT), is capable of retrieving velocity fields in an Area Of Interest (AOI) by using input data ranging from high-resolution slope angle maps derived from Airborne Lidar [22], to medium resolution Landsat-7 or Landsat-8 [23], high resolution satellite images [24,25] and X-band Synthetic Aperture Radar amplitude images [26]. In principle, the minimum sensed displacement corresponds to the spatial resolution of the dataset, but this can be upgraded to subpixel resolution by redundant processing [24] or by pre-processing [27]. With the advent of the Sentinel-2 constellation of multi-sensors satellites with a revisit time of 5 days (that can be even less using the overlap between swaths from adjacent orbits), DIC applications are being developed in order to exploit the potentialities of this continuously increasing dataset [28,29,30]. In general, DIC algorithms are based on two different implementations: (i) the Normalized Cross-Correlation (NCC) algorithm, which operates in the spatial domain [24,31]; (ii) the Fast Fourier Transformation (FFT) approach, which operates in the frequency domain [32,33].
In this paper, we deal with the application of an NCC algorithm to 51 Sentinel-2 Multi Spectral Instrument (MSI) Level-2A acquisitions, to retrieve the spatial and temporal distribution of slope movements in the Ca’ Lita landslide (northern Apennines, Italy) over the period 19 February 2016 to 16 July 2019. Results are integrated with, and validated by, data from in situ continuous GNSS systems, in order to highlight the advantages and limitations of the NCC derived displacement time series for monitoring the spatial distribution of slope movements in moderately rapid landslides.

2. Case Study

The Ca’ Lita landslide is located in the northern Apennines of Italy, in the Emilia-Romagna Region (Figure 1). The landslide is a roto-translational rock slide in the head zone, where flysch rock masses are involved, and evolves a downslope as a translational earthslide-earthflow affecting clayey complexes and debris material from the degradation of flysch rock masses [34]. The total length is approximately 3 km and the maximum width is about 1.4 km, with a maximum depth of the sliding surface of about 42 m measured by inclinometers in the head zone [35].
In April 2004, the Ca’ Lita landslide underwent a first paroxysmal reactivation that caused significant retrogression and advancement of the landslide [36,37]. After that, the landslide was subject to the design and implementation of structural mitigation works [34,35,36,37], consisting of draining trenches, draining wells, check dams and pile-funded anchored retaining walls, which were mostly completed by 2011. The implementation of a monitoring network associated to deep drainage systems evidenced the influence of deep-groundwater inflow in the evolution of the landside [38]. In March 2016, a new paroxysmal reactivation of the landslide took place, resulting in displacements up to tens of meters in the upper landslide and up to more than hundred meters in the earthflow zone. The structural mitigation measures were severely damaged; retaining walls were overturned or compromised and part of the deep-drainage wells lost functionality. From a few days after reactivation, i.e., from 24 March 2016 until June 2016, a GNSS monitoring system, composed of three rover units and one reference master unit [1], was installed in order to support civil protection activities (Figure 2) (rovers of this first monitoring period are indicated with the subscript “a”). The deceleration of the landslide during the following months was recorded by GNSS rovers. By the end of April 2019, cumulated planimetric displacements in the order of 44 m in the upper portion of the earthslide-earthflow (ROV2a) and 2.7 m in the lower portion of the earthslide-earthflow (ROV1a) were reached. In this same period, the roto-translational slide in the head zone did not show relevant movements (ROV3a). In December 2017 to April 2018, more than 40 m of movement took place in the earthslide-earthflow zone. This is known on the basis of the GNSS systems reinstalled in December 2016 that is still operative today (rovers in this second monitoring period are referred to with the subscript “b”). In May to June 2019, another large-scale paroxysmal reactivation took place. In the head zone, 9.6 m of planar and −2 m of vertical displacements were reached in the course of 12 h (ROV3b). Soon after, movements propagated downslope, cumulating up to 70 m planimetric displacement in the upper track zone (ROV2b) and up to 40 m in the lower track zone (ROV1b). The specific characteristics of the GNSS dataset will be presented in details in the “Methods and Data” paragraph.

3. Methods and Data

In this work, the Normalized Cross Correlation procedure by Kääb, Vollmer and Heid [40,41], implemented in the Correlation Image Analysis (CIAS) software has been used. The CIAS software works with pairs of precisely co-registered images in geotiff or tiff-world format acquired from any sensors. Images must have equal vertical and horizontal resolution and need to be single channel 8-bit. The dataset used in this work includes 51 Sentinel-2 Multi Spectral Instrument (MSI) Level-2A acquisitions with a median temporal granularity of 21 days from February 2016 to July 2019 (see Table 1 for details and Figure 3a for processing flowchart). Such data, freely distributed through the SciHub portal of the European Space Agency (ESA) [42], are ortho-images in UTM/WGS84 projection provided as Bottom of Atmosphere (BOA) reflectance images derived from the associated Level-1C products [43]. In this work, data belonging to band 4 of the MSI sensor, characterised by a 10 m spatial resolution, have been used. By iteratively applying the NCC analysis to subsequent pairs of images, a total of 50 displacement maps have been obtained, from which displacement time-series could be derived for every pixel within the AOI. By taking into consideration the dimensions of the landslide (minimum width of 150 m in the track zone) and the expected maximum displacement between two scenes, the settings adopted in CIAS have been a reference block size of 100 m and search area size of 150 m, which allow for maximum theorical planimetric displacements between the two analysed scenes of 70 m. Prior to NCC analysis, the co-registration pre-analysis routine inside CIAS has been run using 14 points outside the landslide as stable benchmarks (Figure 2). More in detail, these benchmarks are represented by elements at ground height (in order to avoid any parallax bias) such as road intersections and corners of cultivation fields located outside the landslide that show no evidence of movements. This allows us to bypass limitations related to the absolute geolocation uncertainty of 20 m without Ground Control Points (GCP) and 12.5 m with GCPs and a 0.3 pixel multi-temporal registration [43].
Continuous GNSS receivers have been operated in the landslide during most of the period analyses with the NCC of Sentinel 2 datasets. GNSS receivers were located at the head zone (ROV3), the upper track zone (ROV2) and the lower track zone (ROV1). The GNSS monitoring dataset covers two time intervals: period “a” from 24 March 2016 to 1 June 2016; period “b” from 28 November 2016 to 15 July 2019 (Figure 3b). During period “a”, the GNSS monitoring array was composed by: an external reference station (double frequency Leica L1-L2 receiver mod. GMX902), and three monitoring rovers (Leica L1 receivers mod. GMX901) indicated as ROV1a, ROV2a and ROV3a. During period “b” the GNSS monitoring array was composed by: two reference stations (a Leica L1-L2 receiver mod. GMX902 in exactly the same position as during period “a”, plus an Emlid L1 receiver model Reach RS installed right next to it, Figure 2a), and three monitoring rovers (two Leica L1 receivers mod. GMX901, indicated as ROV1b and ROV2b, plus an Emlid L1 receiver model “Reach RS” indicated as ROV3b) (Figure 4). With the Leica GNSS array, positioning solutions are computed on the basis of rinex products of different duration (1 h, 30 min and 10 min), and are post-processed by means of the “Leica Spider Software” installed in an industrial PC in the field. Position data are then sent to a remote data server. With the Emlid GNSS array, the ROV3b computes 5 Hz Real Time Kinematic (RTK) solutions by using Radio Technical Commission for Maritime Services (RTCM) stream corrections sent by the Emlid reference station. RTK positioning solutions are calculated by the receiver, and hourly positioning solutions are sent to a remote data server. This type of RTK array, is a cost-effective alternative to the use of high-end Geodetic receivers [44]. A detailed description of the GNSS positioning dataset and of the technical characteristics of the arrays is given in Mulas et al. [39]. The monthly stability check of the reference station is presented in Appendix A.
The analysis of the DIC displacement time-series (DICdispl) has been carried out by comparing them to the planimetric continuous GNSS displacement time-series (GNSSdispl), derived by computing for each rover the median value of GNSS positioning solutions on 1-day long time windows. The DICdispl time series have been retrieved by computing, for each Sentinel-2 acquisition time-step, the mean value of DIC displacement (direction and magnitude) in 60 × 60 m computation matrices centred on the median coordinate values of each GNSS rover. Consequently, the DICdispl time series have been created for the head zone (ROV3a and ROV3b), upper track zone (ROV2a and ROV2b) and lower track zone (ROV1a and ROV1b). The discrepancy between the DICdispl and GNSSdispl time series at variable minimum displacement thresholds has been quantified using the relative % error analysis and the linear fit between time series.

4. Results

By applying NCC to Sentinel 2 data, a total of 50 deformation maps of the Ca’ Lita landslide have been obtained for the period February 2016 to July 2019 with a median periodicity of 21 days (Table 1). The resulting displacements maps allow us to identify changes occurring in the AOI, such as landslide activity or the Secchia river channel migration (Figure 5). We recall that Sentinel-2 scenes were selected by discarding images with cloud cover above the active landslide area. Despite that, cloud cover effect in the remaining frame is still sometimes an issue (Figure 5d). Taking this into consideration, and in order to focus the attention on the Ca’ Lita landslide dynamics, in Figure 6 and Figure 7, we have masked the NCC results outside the active landslide outline.
The displacements map referring to February–March 2016, bracketing the onset of the March 2016 reactivation event, is presented in Figure 6 together with “rose-wind like” plots of displacement vectors, to highlight preferential directions of movement within morphological zones. The map (masked outside the active landslide) evidences higher displacement rates in the head zone and in the higher track zone. Conversely, the lower track zone did not experience relevant displacements. The rose diagrams show that the recorded directions of DIC displacement vectors are coherent with the expected kinematic of the landslide in the different zones. A selection of displacement maps referring to different activity phases in the entire analysis period is presented in Figure 7: the paroxysmal reactivation phase of March 2016, the subsequent suspended movements phase of Summer 2016, the partial acceleration of March 2018 and, finally, the last large-scale reactivation phase of May to June 2019.
The results of the DICdispl time-series extraction with the 60 × 60 m computation matrices centred on the median coordinates values of each GNSS rover are presented in Figure 8. An issue with the construction of the DICdispl time-series is that errors can cumulate at each time step resulting, in the long-term, in unreliable cumulated displacement values. In order to assess this issue, prior to DICdispl time-series extraction, we have filtered the DIC movement maps by using different minimum displacement thresholds ranging from 1 to 10 m. This has allowed us to generate different DICdispl time-series, at variable filtering thresholds, for slope movements around the GNSS (Figure 8). Results show, qualitatively, that the higher the displacement threshold is, the better is the fit between the DICdispl and the GNSSdispl time-series. For instance, by adopting a displacement threshold of ≥10 m (i.e., the pixel size of Sentinel data), the DICdispl over the entire period of analysis reach totals of around 290 m in the upper track zone, 90 m in the lower track zone and 20 m in the head zone, which are quantitatively comparable with GNSSdispl results. On the other hand, by using minimum displacement thresholds ˂10 m (i.e., sub-pixel accuracies), the DICdispl time series tend to overestimate movements, as errors are cumulated over time. In the high mobility GNSS rovers ROV1 and ROV2, the offset between the DICdispl and the GNSSdispl time-series becomes quantitatively significant for minimum displacement thresholds lower than 5 to 7 m. In the low mobility ROV3, the discrepancy is already evident for minimum displacement thresholds lower than 8 m.
A quantitative assessment of the errors associated with the DICdispl time-series retrieval at sub-pixel resolution has been performed by computing the relative % error between DICdispl and GNSSdispl for each minimum displacement threshold considered. The plot in Figure 9a evidences that relative errors are way above 100% for minimum displacement thresholds lower than 2 m (i.e., 0.2 pixels), while they are reduced to less than 15% for minimum displacement thresholds higher than 6 m. In Figure 9b, the coefficients of determination (R2) for the linear fits between DICdispl and GNSSdispl at different minimum displacement thresholds are presented together with the resulting equations (y1 to y10 in relation to thresholds from 1 to 10 m). The R2 indexes are all relatively good, between 0.90 and 0.93, evidencing that even if the relative errors might be very high, the discrepancy between the absolute value of the measured movements is limited even at sub-pixel resolutions. The frequency distribution of relative errors associated with DICdispl at different minimum displacement thresholds is presented in Figure 9c. The results are that the median values of relative errors vary in three groups inside a quite small range. One group is that of results associated to thresholds <3 m, with the higher median relative errors and the larger interquartile ranges (IQR). Another group is associated to thresholds from 3 to 7 m, with slightly lower median error values and a smaller IQR. Finally, a third group is that of results for thresholds higher than 7 m, with the lowest median relative error values and, also, the lowest IQR.

5. Discussion

The set of results obtained indicates that this technique can be a useful tool to monitor the dynamics of large-scale slope movements characterised by moderate velocity ranges. With this respect, the following issues can be discussed.
  • Quantitative errors associated to displacement values retrieved by NCC. Our results prove that NCC movement estimates between two subsequent scenes are: (i) affected by large relative errors for movements lower than 0.3 pixels; (ii) affected by significant errors—yet indicative of ongoing slope dynamics trends—for movements between 0.3 to 0.7 pixels; (iii) quite reliable for movements larger than 0.7 pixels, providing time-series that can be used to monitor slope movements evolution over long periods of time, with median relative errors in the range of 11.66%. The error analysis shows that the CIAS algorithm applied to Sentinel-2 band 4 scenes reaches good performances (approx. 10% error) when displacements are of the same order of the scene geometric resolution (10 m). These results are somehow at odds with the conclusion of other authors, which have reported higher sub-pixel accuracies. Bickel et al. [16], by presenting displacement maps obtained by DIC coherent with ones derived by using GNSS, argue that the velocity of the observed displacement did not seem to have an influence on the correlation accuracy. Mazzanti at al. [25] report sub-pixel accuracies using a DIC methodology based on redundant correlations with multiple master images, and declare precisions up to 1% of the pixel size. Stumpf et al. [30] also reports accuracies in between 1/9th and 1/20th of a pixel by using a two-step method that calibrates 1st grade DICdispl residuals against ground measurements using High Performance Computing (HPC) infrastructures.
  • Applicability of NCC to large displacements and moderate velocity ranges. Results prove that NCC goes beyond the range of landslide velocities that are normally retrieved with the application of other satellite-based techniques, such as Synthetic Aperture Radar (SAR) interferometry. This study indicates that it is possible to detect displacements during activation phases of large landslides. This is true even with displacements of several tens of meters and velocities in ranges from slow to moderate according to landslide movement rate scale by Cruden and Varnes [45]. Similar findings were achieved by other authors in monitoring ice flow measurements in glaciers [28].
  • Identification of spatially distributed displacement patterns. Another advantage is that the output of DIC is spatial displacements maps that can reveal interesting information on the reactivation patterns of large landslides that otherwise are generally monitored in a few specific points only. Results shown in Figure 5 demonstrate that outside the landslide area, the movements computed by the DIC algorithm are substantially nil, strengthening the interpretation of movement patterns retrieved inside the landslide area.
  • Integration with in situ monitoring systems. As a matter of fact, when displacement rates are significant, the maintenance in operation of in situ-systems is quite demanding and potentially very costly (as GNSS or other systems can possibly be lost during the landslide activity). In practice, with NCC, we have been able to obtain long-term time-series that can be used to extend existing—often discontinuous—monitoring time-series from in situ instruments.
  • Main limitations in the application of DIC techniques to Sentinel-2 scenes and the individuation of alternative sensing solutions. One limitation that it is worthwhile to recall is that the MSI onboard the Sentinel-2 constellation is susceptible to cloud cover in the AOI. Even if the 5-day acquisition frequency can help to reduce this issue, this limits the number of MSI acquisitions that can effectively be used in order to retrieve displacements rates by the application of DIC techniques. Even MSI acquisition that is cloud-free over the landslide area can, however, produce false displacement information in areas outside the landside where cloud coverage exists (see for example Figure 5d). Moreover, the presence of snow cover also inhibits the application of DIC algorithms, making the identification of movements during snow-cover or snow melt phases quite problematic. Such issue could be tackled by using satellite constellation with sub metric resolution and up to 12 h revisit time (i.e., “Planet Scope” from “planet” company) [25]. An alternative solution to reduce the influence of cloud cover could be represented by the use of the SAR amplitude scenes. The use of high-resolution COSMO-SkyMed or TerraSAR-X X-band SAR constellations has been already proven by several authors by applying DIC or radargrammetry algorithms [17,26]. Despite the technical advantages of X-band SAR and of high spatial and temporal resolution earth observing satellite constellations, we should recall that they are commercial products while Sentinel-2 scenes are freely distributed by the ESA.
  • Possibility to identify other geomorphic processes using DIC algorithms. The application to the entire AOI has highlighted that the adopted technique can also be used to reveal river dynamics. Results in Figure 5 evidence changes along the Secchia river, which can be ascribed to bedload transport and channel migration. The analysis of such fluvial dynamics goes beyond the scope of this paper.

6. Conclusions

This paper presents a successful application of DIC applied to Sentinel-2 MSI scenes to retrieve slope movement maps and displacement time-series for a large-scale moderate-velocity landslide during a more than three years period. During this time interval, three major reactivations have been back analysed, which resulted in cumulated displacements up to 300 m in the upper track zone. The relative discrepancy between displacements assessed using DIC compared to continuous GNSS measurements has an 11.66% median value using minimum displacement resolutions of 0.7 pixels. This allows us to retrieve time-series with a good fit with displacements obtained by in situ GNSS sensors. In perspective, thanks to the availability of data acquired by Sentinel-2 and Sentinel-1 constellations with high temporal frequency and global coverage, the integrated use of interferometry based (differential interferometry and multi-temporal interferometry) and DIC processing techniques can allow us to bypass flaws of both techniques by fully exploiting the specific ranges of movement measurement capabilities. In the case of the Ca’ Lita landslide, DIC can be applied to monitor rapid accelerations that cannot be resolved by interferometry based techniques, while SAR interferometry can be used to monitor slow movements phases. For future developments, the application of redundant NCC algorithms [25,30] to the Sentinel-2 stack used in this work would be of great interest in order to measure the accuracy improvement that could be achieved by using these techniques. Future research will also focus on the error estimation of DIC applied to Landsat-7 scenes (15 m, panchromatic band) in order to verify the feasibility of analysing historical datasets and following the evolution of large landslides along decades.

Author Contributions

A.C. and M.M. conceived and designed the experiments; A.C., G.C., G.T. and M.M. maintained the in-situ monitoring system; M.M. analysed the data; A.C. supervised the project; M.M. wrote the paper draft; A.C., G.C., G.T. and M.M. revised the manuscript. All authors have read and agreed to the published version of the manuscript.

Funding

This work was supported by the Agency for Civil Protection and Territorial Security of Emilia-Romagna Region, under the framework of the project “Prosecuzione dello studio e del monitoraggio frane Sassi Neri (Farini), in continuo e Pergalla (Bettola) in periodico” (SMART C.I.G.—SIMOG Z5C29E7E0C), activity 4: “Landslide Monitoring with in-situ Permanent GNSS arrays” referring to the permanent monitoring in other sites of interest to the Agency. Principal Investigators: A. Corsini, M. Mulas, G. Ciccarese.

Conflicts of Interest

The authors declare no conflict of interest.

Appendix A

In order to verify the stability of the area where both the GNSS reference stations are located, coordinates of the Leica station have been periodically computed. Differential solutions with respect to a permanent GNSS station belonging to the “Rete Integrata Nazionale GNSS” (RING) national network [46] are obtained. More in detail, daily solutions have been processed, in the period between December 2016 and July 2019, on a monthly basis using as reference the MODE permanent station (30 km NE from Ca’ Lita landslide). Solutions (Figure A1) show a small standard deviation (around 0.01 m for the planimetric components and 0.02 m for the elevation) with respect to their median value (Table A1).
Remotesensing 12 02605 g0a1 550
Figure A1. Stability analysis of the Leica GNSS reference station. Daily coordinate solutions computed monthly from December 2016 to July 2019.
Figure A1. Stability analysis of the Leica GNSS reference station. Daily coordinate solutions computed monthly from December 2016 to July 2019.
Remotesensing 12 02605 g0a1
Table
Table A1. Median value of the coordinate solutions of the Leica reference station.
Table A1. Median value of the coordinate solutions of the Leica reference station.
East (m) 1North (m) 1Elevation (m) 1
628,981.2434,924,121.621499.776
1 WGS84/UTM32N (EPSG: 32632) coordinate system.

References

  1. Corsini, A.; Bonacini, F.; Mulas, M.; Ronchetti, F.; Monni, A.; Pignone, S.; Primerano, S.; Bertolini, G.; Caputo, G.; Truffelli, G.; et al. A portable continuous GPS array used as rapid deployment monitoring system during landslide emergencies in Emilia Romagna. Rend. Online Della Soc. Geol. Ital. 2015, 35, 89–91. [Google Scholar] [CrossRef]
  2. Corsini, A.; Bonacini, F.; Deiana, M.; Giusti, R.; Russo, M.; Ronchetti, F.; Cantini, C. A wireless crackmeters network for the analysis of rock falls at the Pietra di Bismantova natural heritage site (Northern Apennines, Italy). In Proceedings of the Landslides and Engineered Slopes. Experience, Theory and Practice, Napoli, Italy, 12–19 June 2016; pp. 685–690. [Google Scholar]
  3. Mikkelsen, P.E. Field instrumentation. In Landslides Investigation and Mitigation, Transportation Research Board; Turner, A.K., Schuster, R.L., Eds.; Academy Press: Washington, DC, USA, 1996. [Google Scholar]
  4. Mulas, M.; Marnas, M.; Ciccarese, G.; Corsini, A. Sinusoidal wave fit indexing of irreversible displacements for crackmeters monitoring of rockfall areas: Test at Pietra di Bismantova (Northern Apennines, Italy). Landslides 2020, 17, 231–240, Correction to: Sinusoidal wave fit indexing of irreversible displacements for crackmeters monitoring of rockfall areas: Test at Pietra di Bismantova (Northern Apennines, Italy). Landslides 2020, 17, 241–242. [Google Scholar] [CrossRef]
  5. Tonnellier, A.; Helmstetter, A.; Malet, J.P.; Schmittbuhl, J.; Corsini, A.; Joswig, M. Seismic monitoring of soft-rock landslides: The Super-Sauze and Valoria case studies. Geophys. J. Int. 2013. [Google Scholar] [CrossRef] [Green Version]
  6. Barbarella, M.; Fiani, M.; Lugli, A. Landslide monitoring using multitemporal terrestrial laser scanning for ground displacement analysis. Geomat. Nat. Hazards Risk 2015, 6, 398–418. [Google Scholar] [CrossRef]
  7. Corsini, A.; Bonacini, F.; Mulas, M.; Petitta, M.; Ronchetti, F.; Truffelli, G. Long-Term Continuous Monitoring of a Deep-Seated Compound Rock Slide in the Northern Apennines (Italy). In Engineering Geology for Society and Territory-Volume 2 SE-235; Lollino, G., Giordan, D., Crosta, G.B., Corominas, J., Azzam, R., Wasowski, J., Sciarra, N., Eds.; Springer International Publishing: Cham, Switzerland, 2015; pp. 1337–1340. ISBN 978-3-319-09056-6. [Google Scholar]
  8. Hu, J.; Guo, J.; Xu, Y.; Zhou, L.; Zhang, S.; Fan, K. Differential Ground-Based Radar Interferometry for Slope and Civil Structures Monitoring: Two Case Studies of Landslide and Bridge. Remote Sens. 2019, 11, 2887. [Google Scholar] [CrossRef] [Green Version]
  9. Corsini, A.; Mulas, M. Use of ROC curves for early warning of landslide displacement rates in response to precipitation (Piagneto landslide, Northern Apennines, Italy). Landslides 2016, 14, 1241. [Google Scholar] [CrossRef]
  10. Corsini, A.; Castagnetti, C.; Bertacchini, E.; Rivola, R.; Ronchetti, F.; Capra, A. Integrating airborne and multi-temporal long-range terrestrial laser scanning with total station measurements for mapping and monitoring a compound slow moving rock slide. Earth Surf. Process. Landf. 2013, 38, 1330–1338. [Google Scholar] [CrossRef]
  11. Thiebes, B.; Tomelleri, E.; Mejia-Aguilar, A.; Rabanser, M.; Schlögel, R.; Mulas, M.; Alessandro, C. Assessment of the 2006 to 2015 Corvara Landslide Evolution Using a UAV-Derived DSM and Orthophoto. In Proceedings of the Landslides and Engineered Slopes. Experience, Theory and Practice, Napoli, Italy, 12–19 June 2016; pp. 1897–1902. [Google Scholar]
  12. Niethammer, U.; James, M.R.; Rothmund, S.; Travelletti, J.; Joswig, M. UAV-based remote sensing of the Super-Sauze landslide: Evaluation and results. Eng. Geol. 2012, 128, 2–11. [Google Scholar] [CrossRef]
  13. Rossini, M.; Di Mauro, B.; Garzonio, R.; Baccolo, G.; Cavallini, G.; Mattavelli, M.; De Amicis, M.; Colombo, R. Rapid melting dynamics of an alpine glacier with repeated UAV photogrammetry. Geomorphology 2018, 304, 159–172. [Google Scholar] [CrossRef]
  14. Karantanellis, E.; Marinos, V.; Vassilakis, E.; Christaras, B. Object-Based Analysis Using Unmanned Aerial Vehicles (UAVs) for Site-Specific Landslide Assessment. Remote Sens. 2020, 12, 1711. [Google Scholar] [CrossRef]
  15. Valkaniotis, S.; Papathanassiou, G.; Ganas, A. Mapping an earthquake-induced landslide based on UAV imagery; case study of the 2015 Okeanos landslide, Lefkada, Greece. Eng. Geol. 2018, 245, 141–152. [Google Scholar] [CrossRef]
  16. Bickel, V.; Manconi, A.; Amann, F. Quantitative Assessment of Digital Image Correlation Methods to Detect and Monitor Surface Displacements of Large Slope Instabilities. Remote Sens. 2018, 10, 865. [Google Scholar] [CrossRef] [Green Version]
  17. Fratarcangeli, F.; Nascetti, A.; Capaldo, P.; Mazzoni, A.; Crespi, M. Centimeter COSMO-Skymed range measurements for monitoring ground displacements. Int. Arch. Photogramm. Remote Sens. Spat. Inf. Sci. ISPRS Arch. 2016, 41, 815–820. [Google Scholar] [CrossRef]
  18. Mulas, M.; Petitta, M.; Corsini, A.; Schneiderbauer, S.; Mair, V.; Iasio, C. Long-term monitoring of a deep-seated, slow-moving landslide by mean of C-band and X-band advanced interferometric products: The Corvara in Badia case study (Dolomites, Italy). ISPRS Int. Arch. Photogramm. Remote Sens. Spat. Inf. Sci. 2015, XL-7/W3, 827–829. [Google Scholar] [CrossRef] [Green Version]
  19. Schlögel, R.; Thiebes, B.; Mulas, M.; Cuozzo, G.; Notarnicola, C.; Schneiderbauer, S.; Crespi, M.; Mazzoni, A.; Mair, V.; Corsini, A. Multi-Temporal X-Band Radar Interferometry Using Corner Reflectors: Application and Validation at the Corvara Landslide (Dolomites, Italy). Remote Sens. 2017, 9, 739. [Google Scholar] [CrossRef] [Green Version]
  20. Corsini, A.; Farina, P.; Antonello, G. Space-borne and ground-based SAR interferometry as tools for landslide hazard management in civil protection. Int. J. Remote Sens. 2006, 27, 2351–2369. [Google Scholar] [CrossRef]
  21. Bardi, F.; Raspini, F.; Frodella, W.; Lombardi, L.; Nocentini, M.; Gigli, G.; Morelli, S.; Corsini, A.; Casagli, N. Monitoring the rapid-moving reactivation of earth flows by means of GB-InSAR: The April 2013 Capriglio Landslide (Northern Appennines, Italy). Remote Sens. 2017, 9, 165. [Google Scholar] [CrossRef] [Green Version]
  22. Daehne, A.; Corsini, A. Kinematics of active earthflows revealed by digital image correlation and DEM subtraction techniques applied to multi-temporal LiDAR data. Earth Surf. Process. Landf. 2013, 38, 640–654. [Google Scholar] [CrossRef]
  23. Caporossi, P.; Mazzanti, P.; Bozzano, F. Digital Image Correlation (DIC) Analysis of the 3 December 2013 Montescaglioso Landslide (Basilicata, Southern Italy): Results from a Multi-Dataset Investigation. ISPRS Int. J. Geo-Inf. 2018, 7, 372. [Google Scholar] [CrossRef] [Green Version]
  24. Stumpf, A.; Malet, J.P.; Delacourt, C. Correlation of satellite image time-series for the detection and monitoring of slow-moving landslides. Remote Sens. Environ. 2017, 189, 40–55. [Google Scholar] [CrossRef]
  25. Mazzanti, P.; Caporossi, P.; Muzi, R. Sliding Time Master Digital Image Correlation Analyses of CubeSat Images for landslide Monitoring: The Rattlesnake Hills Landslide (USA). Remote Sens. 2020, 12, 592. [Google Scholar] [CrossRef] [Green Version]
  26. Mulas, M.; Corsini, A.; Cuozzo, G.; Callegari, M.; Thiebes, B.; Mair, V. Quantitative monitoring of surface movements on active landslides by multi-temporal, high-resolution X-Band SAR amplitude information: Preliminary results. In Proceedings of the Landslides and Engineered Slopes. Experience, Theory and Practice, Napoli, Italy, 12–19 June 2016; Volume 3, pp. 1511–1516. [Google Scholar]
  27. Debella-Gilo, M.; Kääb, A. Sub-pixel precision image matching for measuring surface displacements on mass movements using normalized cross-correlation. Remote Sens. Environ. 2011, 115, 130–142. [Google Scholar] [CrossRef] [Green Version]
  28. Kääb, A.; Winsvold, S.; Altena, B.; Nuth, C.; Nagler, T.; Wuite, J. Glacier Remote Sensing Using Sentinel-2. Part I: Radiometric and Geometric Performance, and Application to Ice Velocity. Remote Sens. 2016, 8, 598. [Google Scholar] [CrossRef] [Green Version]
  29. Manconi, A.; Kourkouli, P.; Caduff, R.; Strozzi, T.; Loew, S. Monitoring Surface Deformation over a Failing Rock Slope with the ESA Sentinels: Insights from Moosfluh Instability, Swiss Alps. Remote Sens. 2018, 10, 672. [Google Scholar] [CrossRef] [Green Version]
  30. Stumpf, A.; Michéa, D.; Malet, J.P. Improved co-registration of Sentinel-2 and Landsat-8 imagery for Earth surface motion measurements. Remote Sens. 2018, 10, 160. [Google Scholar] [CrossRef] [Green Version]
  31. Kääb, A. Monitoring high-mountain terrain deformation from repeated air- and spaceborne optical data: Examples using digital aerial imagery and ASTER data. ISPRS J. Photogramm. Remote Sens. 2002, 57, 39–52. [Google Scholar] [CrossRef]
  32. Ayoub, F.; Leprince, S.; Avouac, J.P. User’s Guide to COSI-Corr: Co-Registration of Optically Sensed Images and Correlation; California Institute of Technology: Pasadena, CA, USA, 2017. [Google Scholar]
  33. Guizar-Sicairos, M.; Thurman, S.T.; Fienup, J.R. Efficient subpixel image registration algorithms. Opt. Lett. 2008, 33, 156. [Google Scholar] [CrossRef] [Green Version]
  34. Borgatti, L.; Corsini, A.; Barbieri, M.; Sartini, G.; Truffelli, G.; Caputo, G.; Puglisi, C. Large reactivated landslides in weak rock masses: A case study from the Northern Apennines (Italy). Landslides 2006, 3, 115–124. [Google Scholar] [CrossRef]
  35. Corsini, A.; Borgatti, L.; Caputo, G.; De Simone, N.; Sartini, G.; Truffelli, G. Investigation and monitoring in support of the structural mitigation of large slow moving landslides: An example from Ca’Lita (Northern Apennines, Reggio Emilia, Italy). Nat. Hazard Earth Sci. 2006, 6, 55–61. [Google Scholar] [CrossRef]
  36. Borgatti, L.; Corsini, A.; Marcato, G.; Ronchetti, F.; Zabuski, L. Appraise the structural mitigation of landslide risk via numerical modelling: A case study from the northern Apennines (Italy). Georisk Assess. Manag. Risk Eng. Syst. Geohazards 2008, 2, 141–160. [Google Scholar] [CrossRef]
  37. Cervi, F.; Ronchetti, F.; Martinelli, G.; Bogaard, T.A.; Corsini, A. Origin and assessment of deep groundwater inflow in the Ca’ Lita landslide using hydrochemistry and in situ monitoring. Hydrol. Earth Syst. Sci. 2012, 16, 4205–4221. [Google Scholar] [CrossRef] [Green Version]
  38. Ronchetti, F.; Borgatti, L.; Cervi, F.; Corsini, A. Hydro-mechanical features of landslide reactivation in weak clayey rock masses. Bull. Eng. Geol. Environ. 2010, 69, 267–274. [Google Scholar] [CrossRef]
  39. Mulas, M.; Ciccarese, G.; Truffelli, G.; Corsini, A. Displacements of an active moderately rapid landslide–A dataset retrieved by continuous GNSS arrays. Data 2020, 5, 71. [Google Scholar] [CrossRef]
  40. Kääb, A.; Vollmer, M. Surface geometry, thickness changes and flow fields on creeping mountain permafrost: Automatic extraction by digital image analysis. Permafr. Periglac. Process. 2000, 11, 315–326. [Google Scholar] [CrossRef] [Green Version]
  41. Heid, T.; Kääb, A. Evaluation of existing image matching methods for deriving glacier surface displacements globally from optical satellite imagery. Remote Sens. Environ. 2012, 118, 339–355. [Google Scholar] [CrossRef]
  42. ESA Copernicus Open Access Hub. Available online: https://scihub.copernicus.eu/dhus/#/home (accessed on 23 March 2020).
  43. ESA Sentinel-2 MSI User Guide. Available online: https://earth.esa.int/web/sentinel/user-guides/sentinel-2-msi (accessed on 23 March 2020).
  44. Madson, A.; Fielding, E.; Sheng, Y.; Cavanaugh, K. High-Resolution Spaceborne, Airborne and in Situ Landslide Kinematic Measurements of the Slumgullion Landslide in Southwest Colorado. Remote Sens. 2019, 11, 265. [Google Scholar] [CrossRef] [Green Version]
  45. Cruden, D.N.; Varnes, D.J. Landslides types and processes. In Landslides: Investigation and Mitigation; Turner, A.K., Shuster, R.L., Eds.; Transportation Research Board: Washington, DC, USA, 1996; pp. 36–75. [Google Scholar]
  46. Istituto Nazionale di Geofisica e Vulcanologia (INGV). INGV RING Working Group (2016), Rete Integrata Nazionale GPS (RING). Available online: http://ring.gm.ingv.it (accessed on 10 August 2020). [CrossRef]
Remotesensing 12 02605 g001 550
Figure 1. Ca’ Lita landslide location. Grid coordinates: WGS84/UTM32N (EPSG: 32632).
Figure 1. Ca’ Lita landslide location. Grid coordinates: WGS84/UTM32N (EPSG: 32632).
Remotesensing 12 02605 g001
Remotesensing 12 02605 g002 550
Figure 2. Global Navigation Satellite System (GNSS) records. (a) Planimetric view; Plots (b,d,f) present relative coordinate variations recorded in the period from 24 March 2016 to 15 July 2016. Plots (c,e,g) present relative coordinate variations recorded in the period from 28 November 2016 to 15 July 2019. GNSS time-series are presented in Mulas et al. [39]. Grid coordinates: WGS84/UTM32N (EPSG: 32632).
Figure 2. Global Navigation Satellite System (GNSS) records. (a) Planimetric view; Plots (b,d,f) present relative coordinate variations recorded in the period from 24 March 2016 to 15 July 2016. Plots (c,e,g) present relative coordinate variations recorded in the period from 28 November 2016 to 15 July 2019. GNSS time-series are presented in Mulas et al. [39]. Grid coordinates: WGS84/UTM32N (EPSG: 32632).
Remotesensing 12 02605 g002
Remotesensing 12 02605 g003 550
Figure 3. (a) Methodological workflow; (b) Temporal-coverage of Sentinel-2 and GNSS data in Ca’ Lita landslide.
Figure 3. (a) Methodological workflow; (b) Temporal-coverage of Sentinel-2 and GNSS data in Ca’ Lita landslide.
Remotesensing 12 02605 g003
Remotesensing 12 02605 g004 550
Figure 4. GNSS receivers installed at the Ca’ Lita landslide. (a) Reference stations. Within the red squares are the Leica post-processing unit (left) and the Leica double frequency receiver on the right (model GMX902). In the blue square is the Emlid single frequency low cost receiver (model Reach RS) acting as a reference station. (b) Reach RS receiver is installed as ROV3b rover unit and its photovoltaic power unit. (c) Leica GMX901 receiver is installed as ROV2b rover unit and its photovoltaic power unit. During period “a”, ROV1a, ROV2a and ROV3a were equipped with Leica GMX901 receivers. During period “b”, ROV1b and ROV2b were equipped with Leica GMX901 receivers, while ROV3b was equipped with the low-cost Emlid Reach RS.
Figure 4. GNSS receivers installed at the Ca’ Lita landslide. (a) Reference stations. Within the red squares are the Leica post-processing unit (left) and the Leica double frequency receiver on the right (model GMX902). In the blue square is the Emlid single frequency low cost receiver (model Reach RS) acting as a reference station. (b) Reach RS receiver is installed as ROV3b rover unit and its photovoltaic power unit. (c) Leica GMX901 receiver is installed as ROV2b rover unit and its photovoltaic power unit. During period “a”, ROV1a, ROV2a and ROV3a were equipped with Leica GMX901 receivers. During period “b”, ROV1b and ROV2b were equipped with Leica GMX901 receivers, while ROV3b was equipped with the low-cost Emlid Reach RS.
Remotesensing 12 02605 g004
Remotesensing 12 02605 g005 550
Figure 5. Changes occurring in the area of interest evidenced by digital image correlation results. (a,b) Landslide activity and Secchia river dynamics. (c) No landslide and river activity. (d) Cloud cover effect. Grid coordinates: WGS84/UTM32N (EPSG: 32632).
Figure 5. Changes occurring in the area of interest evidenced by digital image correlation results. (a,b) Landslide activity and Secchia river dynamics. (c) No landslide and river activity. (d) Cloud cover effect. Grid coordinates: WGS84/UTM32N (EPSG: 32632).
Remotesensing 12 02605 g005
Remotesensing 12 02605 g006 550
Figure 6. Displacement obtained by means of the DIC analysis between Sentinel-2 scenes acquired on 19/02/2016 and 20/03/2016. (a) DICdispl modules within the active landslide area. (1) Head zone; (2) Upper track zone; (3) Lower track zone. (b) Displacement module and geographic direction for all the active landslides (up) and according to the respective geomorphic zone (down). Grid coordinates: WGS84 / UTM32N (EPSG: 32632).
Figure 6. Displacement obtained by means of the DIC analysis between Sentinel-2 scenes acquired on 19/02/2016 and 20/03/2016. (a) DICdispl modules within the active landslide area. (1) Head zone; (2) Upper track zone; (3) Lower track zone. (b) Displacement module and geographic direction for all the active landslides (up) and according to the respective geomorphic zone (down). Grid coordinates: WGS84 / UTM32N (EPSG: 32632).
Remotesensing 12 02605 g006
Remotesensing 12 02605 g007 550
Figure 7. Activity stages of the Ca’ Lita landslide. (a) March 2016 reactivation; (b) Summer 2016 suspended phase; (c) March 2018 reactivation and (d) Spring 2019 reactivation. (1) Head zone. (2) Upper track zone. (3) Lower track zone. Grid coordinates: WGS84/UTM32N (EPSG: 32632).
Figure 7. Activity stages of the Ca’ Lita landslide. (a) March 2016 reactivation; (b) Summer 2016 suspended phase; (c) March 2018 reactivation and (d) Spring 2019 reactivation. (1) Head zone. (2) Upper track zone. (3) Lower track zone. Grid coordinates: WGS84/UTM32N (EPSG: 32632).
Remotesensing 12 02605 g007
Remotesensing 12 02605 g008 550
Figure 8. Comparison between cumulative planimetric DICdispl and GNSSdispl. GNSSdispl time-series of period “b” is translated on the y-axis in order to match the minimum DICdispl value recorded on the first date of the “b” period.
Figure 8. Comparison between cumulative planimetric DICdispl and GNSSdispl. GNSSdispl time-series of period “b” is translated on the y-axis in order to match the minimum DICdispl value recorded on the first date of the “b” period.
Remotesensing 12 02605 g008
Remotesensing 12 02605 g009 550
Figure 9. Error analysis: (a) relative error of DICdispl as a function of GNSSdispl, where the grey area is the 95% confidence band for the red dashed regression curve; (b) DICdispl vs. GNSSdispl: y1 to y10 are the equations resulting from linear fit of points belonging to thresholds from 1 to 10 m and the respective coefficient of determination. Displacements smaller than 1 m are excluded; (c) Percentage error according to displacement threshold classes. Boxplots ranges: extremes of the boxes are the 1st and 3rd quartile (Q1 and Q2). The black lines depict median values, the lower whisker is “Q1-1.5 × IQR” and the upper whisker is “Q2-1.5 × IQR” (where IQR stands for “interquartile range” and is the difference between Q2 and Q1).
Figure 9. Error analysis: (a) relative error of DICdispl as a function of GNSSdispl, where the grey area is the 95% confidence band for the red dashed regression curve; (b) DICdispl vs. GNSSdispl: y1 to y10 are the equations resulting from linear fit of points belonging to thresholds from 1 to 10 m and the respective coefficient of determination. Displacements smaller than 1 m are excluded; (c) Percentage error according to displacement threshold classes. Boxplots ranges: extremes of the boxes are the 1st and 3rd quartile (Q1 and Q2). The black lines depict median values, the lower whisker is “Q1-1.5 × IQR” and the upper whisker is “Q2-1.5 × IQR” (where IQR stands for “interquartile range” and is the difference between Q2 and Q1).
Remotesensing 12 02605 g009
Table
Table 1. List of Sentinel-2 acquisitions.
Table 1. List of Sentinel-2 acquisitions.
YearMonth and Day of AcquisitionTotal
201619 February; 20 March; 19 April; 22 May; 18 June; 28 July; 17 August; 16 October; 8 November; 8 December10
20177 January; 27 January; 26 February; March; 14 April; 14 May; 3 June; 23 June; 13 July; 2 August; 22 August; 21 September; 16 October; 3 November; 20 November; 13 December; 25 December17
201817 January; 23 March; 19 April; 19 May; 28 June; 18 July; 31 July; 20 August; 9 September; 29 September; 24 October; 8 November; 28 November; 8 December14
201914 January; 13 February; 15 March; 30 March; 17 April; 2 May; 3 June; 16 June; 1 July; 16 July10
51

Share and Cite

MDPI and ACS Style

Mulas, M.; Ciccarese, G.; Truffelli, G.; Corsini, A. Integration of Digital Image Correlation of Sentinel-2 Data and Continuous GNSS for Long-Term Slope Movements Monitoring in Moderately Rapid Landslides. Remote Sens. 2020, 12, 2605. https://doi.org/10.3390/rs12162605

AMA Style

Mulas M, Ciccarese G, Truffelli G, Corsini A. Integration of Digital Image Correlation of Sentinel-2 Data and Continuous GNSS for Long-Term Slope Movements Monitoring in Moderately Rapid Landslides. Remote Sensing. 2020; 12(16):2605. https://doi.org/10.3390/rs12162605

Chicago/Turabian Style

Mulas, Marco, Giuseppe Ciccarese, Giovanni Truffelli, and Alessandro Corsini. 2020. "Integration of Digital Image Correlation of Sentinel-2 Data and Continuous GNSS for Long-Term Slope Movements Monitoring in Moderately Rapid Landslides" Remote Sensing 12, no. 16: 2605. https://doi.org/10.3390/rs12162605

APA Style

Mulas, M., Ciccarese, G., Truffelli, G., & Corsini, A. (2020). Integration of Digital Image Correlation of Sentinel-2 Data and Continuous GNSS for Long-Term Slope Movements Monitoring in Moderately Rapid Landslides. Remote Sensing, 12(16), 2605. https://doi.org/10.3390/rs12162605

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop