Advantages of Using UAV Digital Photogrammetry in the Study of Slow-Moving Coastal Landslides

Abstract

The aim of this study is to demonstrate the advantages of using micro drones in the study of large slow-moving landslides, which are widespread along the northwestern coast of Malta. In particular, attention was given to the inventory and analysis of gravity-induced joints and megaclast deposits at four study sites selected due to the presence of remarkable examples of lateral spreads evolving into block slides. The research was carried out by means of Google Earth (GE) image analysis and uncrewed aerial vehicle digital photogrammetry (UAV-DP). UAV-DP outputs enabled the identification and characterization of tens of persistent joints (locally exceeding 150 m) and permitted the size categorization of thousands of blocks. With reference to gravity-induced joints, a favorable agreement was found between existing datasets (mainly based on the integration of GE analysis and field survey) and UAV-DP outputs in terms of the identification of joints and their persistence. Conversely, the use of the UAV-DP technique showed significant advantages in terms of joint aperture determination (even exceeding 1 m) and distribution setting. Regarding the extensive megaclast deposits, UAV-DP enabled the identification of 8943 individuals which, compared with the 5059 individuals identified by GE analysis, showed an increase in the total population of 76%. This is related to the high accuracy of DP-derived orthomosaics and 3D models, which are particularly useful for identifying detached blocks. The inexpensive technique used in this research highlights its potential for being extended to other rocky coastal areas affected by slow-moving landslides.

1. Introduction

Landslides are one of the most frequent geohazards in Europe and can severely impact coastal areas [1,2,3,4,5,6,7,8]. In particular, landslides are common along the European coasts of the Mediterranean Sea and can include different types of transport mechanisms such as falls, slides, spreads, or flows [9,10,11,12,13,14]. Among these, lateral spreading and block sliding processes are widely represented in the Mediterranean coasts, such as in northern Mallorca (Spain) [15], in the western part of Liguria (Northern Italy) [16], Tropea area (Southern Italy) [17], Scopello in northern Sicily [18] and in the Maltese islands [14,19,20].

The term lateral spreading, and more specifically rock spreading, defines a lateral extension in cracked rock masses that can take place in fragile rocks capping visco-plastic terrains [21]. The loading pressure of resistant materials favors the lateral extension of ductile terrains, which are laterally unconfined. This lateral motion of soft terrains generates large bulging processes that form extensive ridges and cause large cracking processes of overlain resistant rock masses. The resulting joints are crucial for the identification of rock spreads. These discontinuities often exceed 100 m in length and are frequently parallelly oriented to the cliff. The fractures can form a broad network of discontinuities and isolate wide sectors of rock masses prone to be detached downslope [22]. For this reason, landslide mapping in areas affected by lateral spreading must include a detailed inventory of gravity-induced joints that are a key factor for the early recognition of this type of mass movement, characterized by a very slow kinematic.

According to Soldati [23] and Pánek and Klimeš [24], lateral spreading can be ascribed to deep-seated gravitational slope deformation (DGSD). DGSDs deform and displace rock masses in extensional, shearing, and mixed modes at very low speeds over long periods [24]. Subsurface geometry is often unknown and DGSDs are characterized by the absence of a continuous or well-defined sliding surface [23]. As a result of these deformations, typical landforms are open tension cracks, grabens and karst-like depressions in brittle rocks and bulges in the underlain material [23,24,25].

Block slide is a slow-moving failure that often affects large coastal areas [1]. Block slides are characterized by markedly noncircular slip surfaces forming a combination of a steep part at the crown and a planar part along the toe area [26]. The slip surface is contained in clays and can also involve a thin limestone caprock. Examples of large block slides with this geological superimposition occur along the Bulgarian coast, where 30–40 m of limestone caps 100 m of clays [27].

In this study, we investigated rock spreads and block slides that severely affect the northwestern coast of Malta using the data collected by a commercial and inexpensive drone. Drones are also known as remotely piloted aircraft systems (RPAS) or UAVs and their use in geosciences has been spreading rapidly, mainly in landslide mapping and geomorphological investigations [28,29]. The use of UAVs has expanded progressively in the last decade related to the fact that commercial drones represent a cheap and rapid solution for the acquisition of detailed images of an area of interest [30].

Typical outputs of drone investigations are high or very high-resolution imagery (HR/VHR), digital elevation models (DEMs), and 3D models [31]. These outcomes are fundamental for desk activities after the completion of field surveys, as they provide support in the identification and mapping of landforms. The generation of DEMs can also be achieved through the use of light detection and ranging (LiDAR) surveys, which have become a fairly common and consolidated procedure in recent years [32]. Typically, there are two main configurations for LiDAR surveys. We refer to airborne laser scanning (ALS) if the sensor is mounted over flying platforms, such as helicopters or airplanes, and terrestrial laser scanning (TLS) if the instrument is steadily positioned over a topographic tripod. ALS surveys can produce highly accurate terrain models over large and remote areas but are extremely expensive and require long and careful planning. In contrast, TLS is much more affordable from an economic viewpoint and better suited for detailed surveys, as the logistics of field measurement campaigns can be very complicated and time-consuming. DEMs derived from LiDAR surveys are certainly more accurate than those generated by the UAV-DP technique, and in landslide investigations, they are mainly used in slope-failure analysis for investigations of steep slopes affected by rock falls or topples [32]. Conversely, UAV-DP is a robust and reliable technique for the investigation of extensive landslides characterized by large transport areas and in all the cases that require a panoramic view from an aerial survey [33]. According to Yao et al. [34], DEMs acquired using drone surveys can provide a centimeter accuracy and for this reason aerial datasets can be used for the investigations of landslide kinematics [35,36,37] and sea cliff stability analysis [38].

In this study, the UAV-DP technique was used to investigate gravity-induced joints and vast block deposits generated by rock spreads and block slides along the Maltese coast, respectively. These joints, formed along the external portions of limestone plateaus by sub-horizontal or very moderate sloping surfaces, were identified, mapped, and classified according to their persistence and aperture using orthomosaics and 3D models produced by the UAV-DP technique.

The spatial distribution and the above-mentioned joint properties were compared with existing datasets derived from recent investigations [39] or with accurate joint data included in the 1:7,500 geomorphological map of the northwestern coast of Malta [40], which, given its extremely high level of detail can be considered a landslide inventory map (LIM).

Finally, we identified and mapped thousands of coastal blocks performing firstly an aerial image interpretation analysis using GE. Most of these blocks were transported by block slides and were categorized using a size classification. These results were compared to outputs made using 3D models and orthomosaics produced by the UAV-DP method in terms of population and size, and successively discussed.

2. Slow-Moving Landslides at the NW Malta Coast

Coastal landslides are widespread along the northern part of the island of Malta (Figure 1) and represent a remarkable geohazard [41,42,43].

A WSW-ENE oriented fault called the Great Fault crosses Malta and can be considered its main structural feature [44]. The area north of the Great Fault is characterized by the presence of structural high and flat areas. This is caused by a horst and graben system, which produces an alternation of WSW-ENE oriented limestone plateaus and narrow valleys [44,45,46].

There are four upland plateaus made up of resistant Upper Coralline Limestone (UCL) Formation [47,48]: (i) Marfa Ridge plateau; (ii) Mellieha ridge; (iii) Badja ridge; and (iv) Ras Il-Pelegrin promontory (Figure 1C). Conversely, the valleys and coastal slopes surrounding limestone plateaus are made up of Blue Clays (BC). BC and UCL occur mainly in the northern sector of the island and are dominant in the northwestern area compared to the older geological rock units, i.e., the Globigerina Limestone (GL) Formation and the oldest Lower Coralline Limestone (LCL) Formation [47,48].

The peculiar geological and structural setting of the northwestern coast of Malta favors the occurrence of abundant and recurrent landslides [14,39,40,49]. Figure 2 shows the landslide distribution and includes a pie chart of different types of mass movements according to the classification of Cruden and Varnes [50].

The presence of BC capped by rigid UCL is a key factor for the development of lateral spreads affecting the external sectors of Marfa Ridge plateau, Anchor Bay, Ras Il-Wahx and Il-Qarraba peninsula (Figure 3A). Lateral spreading is locally controlled by the network of tectonic discontinuities that affect the UCL plateaus [40]. Block slides are also widespread and transport thousands of blocks from the limestone plateau towards the sea (Figure 3B). These vast block deposits extend well below sea level, as reported in recent studies [51,52,53], which constrain their onset to periods when the sea level was much lower than today (up to 130 m), having the landslide deposits partially been drowned as a consequence of post-glacial sea level rise [54].

According to the methodology developed by Bruno and Ruban [55], these slope-failure accumulations can be classified as coastal megaclast deposits, with the limestone blocks normally larger than 1 m in size. Depending on the main axis dimension, this classification subdivides megaclasts into: (i) blocks (1–10 m); (ii) megablocks (10–100 m); and (iii) superblocks when one length exceeds 100 m [55].

Block slides and rock spreads are contained in the category of slow-moving landslides related to the fact that their kinematics are ‘very slow’ or ‘extremely slow’. These slow-moving landslides dominate with respect to other more rapid types of mass movements that can be considered collateral landslides, according to Pasuto and Soldati [56]. Rock falls and rock topples detached from limestone plateaus are strictly connected to the presence of persistent joints originating from rock spreads, whereas earth flows and earth slides are shallow landslides occurring along the clayey slopes and originating mainly from the loading pressure of huge megaclasts moved by block slides, which cause the squeezing of the underlying BC terrains (Figure 3C). A geotechnical description of earth slides and earth flows occurring near Il-Qarraba peninsula was provided by Dykes [57].

3. Study Sites

The chosen sites for drone surveys are: (i) Marfa Ridge; (ii) Anchor Bay; (iii) Ras Il-Wahx headland; and (iv) Il-Qarraba peninsula (Figure 2 and Figure 4).

The study sites were selected based on the following criteria:

• (1). areas characterized by remarkable examples of rock spreads and block slides;

• (2). the presence of gravity-induced joints and/or well-developed network of discontinuities;

• (3). the occurrence of wide population of blocks in the coastal megaclast deposits;

• (4). the availability of published landslide-related datasets;

• (5). guaranteed safety conditions during UAV take-off and landing phases;

• (6). adequate visibility for the pilot during drone flights.

The slow-moving landslides affecting Il-Qarraba and the northern sector of Anchor Bay were proven to be active by long-term monitoring outputs using global navigation satellite system (GNSS) benchmarks operating since 2005 and 2006, respectively [58,59,60]. Evidence of activity of the block slides of Marfa Ridge and Ras-Il-Wahx were detected through advanced SAR interferometric analysis (A-DInSAR) performed using European space missions ERS and ENVISAT datasets (1992–2009) [61]. Moreover, a landslide susceptibility model published by Piacentini et al. [62] classified the entire northern sector of Anchor Bay and large sectors of Ras Il-Wahx and Il-Qarraba peninsula as ′very high susceptibility’ areas. Conversely, Marfa Ridge study area ranges between ′low′ and ′medium′ susceptibility.

3.1. Marfa Ridge

Marfa Ridge is a plateau reaching an altitude of 125 m. It is characterized by a tabular shape limited by 29 m high cliffs. Lateral spreads and block slides mainly affect the western sector, where a well-developed network of joints can be observed on the edge of the limestone plateau [39]. Downslope approximately a thousand of blocks covers the BC slopes, mainly scattered in proximity of this complex discontinuity network (Figure 4A). Conversely, a 132 m long and particularly deep joint is located on the southern part of the limestone plateau, partially filled by fine debris [39]. Its aperture is quite uniform and reflects an initial stage of lateral spreading, proved by the absence of lowering of the external slice of the limestone rock mass [39]. This is due to the occurrence of a 1 km shore platform of GL [39,40,63]. This GL platform avoids the marine erosion of BC slopes and drastically reduces sea water infiltration into BC slopes, slowing down the rate of movement of landslides [61].

3.2. Anchor Bay

From a geomorphological viewpoint, the semi-oval shape of Anchor Bay is related to a structural control due to the presence of two E-W oriented faults [40,60]. At sea level, these faults originate steep UCL cliffs on the southern side and BC slopes in the northern part of Anchor Bay. The above-mentioned BC terrains are backed by an E-W oriented 18 m high escarpment limiting UCL plateau.

In the unconfined inner northern flank of Anchor Bay, also known as Il-Prajjet, the UCL plateau displays typical landforms related to lateral spreading and block sliding (Figure 5). The limestone plateau highlights four main persistent discontinuities. The most persistent joint reaches a length of 95 m [39]. Figure 5 indicates this 2 m open discontinuity, partially filled by debris. This joint is surely generated by lateral spreading and can be classified as a tension crack. The latter partially isolates and lowered a 43,500 m³ limestone slice outlined with red color in Figure 5.

The combination of lateral extension and lowering was clearly inferred by long-term measurements made by GNSS surveys carried out on Benchmark #102 and Benchmark #103 using the fast static (FS) technique over a 4.5-year period [59,60].

At the bottom five megablocks stand in horizontal positions, probably due to the presence of a horizontal slip surface of the block slide positioned inside the BC terrains. Conversely, Megablock #VII is the only megaclast showing a tilled position due to BC squeezing caused by block sliding [39].

Samples for cosmogenic exposure dating collected by the actual escarpment dated back this block slide at ca. 21.7 ka [52], when the sea level was 130 m lower than at present [54].

3.3. Ras Il-Wahx

The structural upland of Badja Ridge terminates westerly at Ras Il-Wahx, which is a headland with high natural interest [64,65,66]. The large block slides surrounding the UCL plateau are active and record the highest deformation rates of the northwestern coast of Malta, according to the A-DInSAR analysis results reported in previous studies [61,62].

These coastal megaclast deposits continue below sea level up to a distance of 700 m from Ras Il-Wahx coastline [51]. Four superblocks are located at the head of block sliding area and are surrounded by tens of megablocks. The uphill-facing position of these superblocks suggests that there is a rotational mechanism of the block slide (Figure 4C). These superblock rotations correspond to multiple rotational slides controlled by deep-seated failure surfaces rooted in BC terrains [39].

Rock spreading involves a narrow sector of Badja Ridge plateau, as indicated by two joints located on the edge of its western part. Their lengths are 50 m and 22 m [39].

3.4. Il-Qarraba Peninsula

The popular inlets of Ghajn Tuffieha Bay and Gnejna Bay are separated by the head-shape peninsula of Il-Qarraba and narrow isthmus that connects it to the mainland (Figure 4D). Il-Qarraba peninsula represents a unique landform of the Maltese archipelago. This is related to the fact that BC terrains are capped in the central part with an 8300 m² hard and jointed UCL caprock. Steep high cliffs limit this tabular caprock, which reaches an altitude of 43 m (Figure 4D and Figure 6) [39].

The western portion of this caprock shows the most extensive and conspicuous evidence of recent deformations caused by slow-moving landslides [39,58]. A well-developed network of joints located on the western side partially isolates five huge rock masses, ranging between 550 m³ and 3300 m³ [39]. Rock mass #II and Rock mass #V are lowered with respect to the caprock, showing the evolution of lateral spreading into block sliding [39]. This network of open, deep, and vertical joints witnesses the occurrence of lateral spreading and suggests an advanced stage of its geomorphological evolution.

Over a thousand of caprock-derived megaclasts surrounds the peninsula, mainly grouped in extensive deposits located in its western and northern sectors (Figure 6). A 35 m megablock and two others smaller megablocks facing the western portion of the caprock cliff indicated that the block slide is a compound slide. This above-cited block is in an uphill-facing position and witnesses a rotational movement connected to the presence of a slip surface. The latter causes rotational deformations of blocks in the inner part of the peninsula and translational displacement in areas close to the sea.

4. Materials and Methods

The investigation of the study sites was carried out using GE images and UAV-DP. GE image analyses enabling megaclast inventory and categorization according to their size were considered functional to subsequent drone surveys. The results of the analysis of already existing datasets and the outputs of GE image analyses were compared with the outcomes of geomorphological investigations using UAV-DP products such as 3D models and orthomosaics for the identification and categorization of gravity-induced joints and coastal megaclast deposits, respectively. The existing analyzed datasets included the inventory of gravity-induced joints and accurate quantitative data of aperture and persistence of each joint [39,40]. The workflow, including a simplified scheme of UAV-DP procedures applied for this study, is summarized in Figure 7.

A micro drone as the DJI Mavic drone^TM (DJI, Nanchan District, Shenzhen, China) was preferred as a good compromise between image quality and flight autonomy, and is capable of automatic planning. One of the vital merits of DJI Mavic drone^TM, when compared to other UAVs, is its small size and light weight, which allows for investigation with a flight time of 30 min and an area of 100,000 m² at a ground resolution of 2 cm/pixel with one battery [67]. Table 1 lists the main technical and operational characteristics of the DJI Mavic^TM.

UAVs are components of an uncrewed aircraft system (UAS), which also includes a remote controller, a tablet for uploading flight plans, and a cable to connect them. This drone is categorized as a very light multicopter (VL/MC) or a micro drone, according to the weight classification proposed by Brooke-Holland [68]. Table 2 lists the drone types of the abovementioned classification.

The DJI Mavic^TM is a vertical takeoff and landing (VTOL) quadricopter [69] and can fly for up to 30 min in idyllic weather conditions. This drone can fly in every direction—vertically or horizontally—and can maintain a fixed position. Furthermore, it supports flight plans computed by software and can navigate automatically. The above-mentioned characteristics make it perfect for landslide areas, which are frequently hard-to-reach zones [70].

4.1. UAV Flight Planning and Image Aquisition

The flight paths were planned using DJI Flightplanner software^TM (Aeroscientific, Blackwood, Australia) and Litchi software^TM (VC Technology, London, UK), on a continuous route at an average speed optimized to guarantee an overlap of 80% on both sizes of the CMOS sensor taking pictures every 2 s without interrupting the flight path. Weather and light condition should be considered in advance to maintain at least a 1/250 s shutter speed at a low ISO setting to avoid motion blur on the pictures. A three-day drone survey was carried out in February 2020.

For safety reasons, the drone flights were planned with particular attention paid to cliffs or any significant variation in altitude. A constant altitude above takeoff elevation was set up over the limestone plateaus using digital terrain models by Google Satellite and uploaded by Litchi software^TM during flight operations, whereas the altitude was higher over block sliding areas. The average ground resolution above downslope terrains is of 3 mm/pixel. These differences exist because a few uphill-facing blocks are rotated from toppling processes reaching a higher altitude than the plateau cliff edge. The DJI Mavic^TM is equipped with sensors and navigated automatically following flight plans computed using DJI Flightplanner software^TM. To improve the matching ability of the software avoiding shooting raw format files, the in-camera shooting parameters were set to the lowest contrast [71]. A total of 2820 photos were taken in separate flights as shown in Table 3.

The quality of the photos was checked during the flight phase and after landing to integrate areas with low-quality pictures. We relied on software capability to detect out-of-focus pictures and remove them.

4.2. 3D Model Production

Postproduction was performed using Agisoft Metashape (Agisoft LLC, St. Petersburg, Russia, 2019) following six steps: (i) image alignment; (ii) point cloud thickening; (iii) mesh generation; (iv) DEM production and exportation; (v) texturing; and (vi) production and exportation of orthomosaics. To analyze the geomorphological features of slow-moving landslides and compare them with existing datasets and GE mapping of megaclasts, DEMs and orthomosaics were imported in QGIS version 3.10 A Coruña (QGIS Development Team, QGIS Geographic Information System, Open Source Geospatial Foundation Project).

To generate a single project file that allows to perform batch processing in several areas, each study site was uploaded as a separate chunk on the Agisoft Metashape^TM software.

Ten different chunks were created; some were joined after the first alignment process. The Marfa Ridge model was produced via the alignment of three different chunks; the first for the UCL plateau, a second for the downslope areas, and a third using pictures taken piloting manually the drone, for assessing the vertical surface [72]. The Il-Qarraba model, on the other hand, was generated in a single chunk to avoid alignment problems due to the low overlap between the two series of images.

4.2.1. Generation of Point Clouds

First, the images were analyzed using the Metashape^TM Estimate Image Quality tool to obtain an average quality value for every single picture [73]. The main goal of this procedure is to detect out-of-focus and negatively exposed pictures. MetashapeTM software permits a selection based on the quality of the pictures. For the first phase, only pictures with a quality marked as lower than 0.7 were removed even if some slightly out-of-focus picture remained in the selection. After the classification and removal of motion-blurred pictures (approximately 0.5%), the alignment task was performed on the entire workspace. Reference preselection enables the use of GPS data provided by the drone to select the image to match; dealing with medium and large numbers of images, as in this work, reference preselection can lower up to 75% of the total time. Accuracy was set to High and no cap to the Key and Tie point numbers were set, allowing the software to detect their maximum amount. Adaptive Camera Model Fitting was enabled, as recommended when working with relatively low-quality non-photographic lenses. Generally, zenithal pictures generate lower quality results in comparison with other photogrammetric acquisitions suitable for smaller objects [74]. Thus, a manual refinement of the cloud was performed to clean the dataset. Using the Gradual Selection tool, approximately 1/3 of the paramount points were removed according to reprojection errors, reconstruction uncertainty and projection accuracy criteria. After this cleaning procedure, the remaining points represent those generated with better confidence [75], with a remaining reprojection error limited, at least to 0.7 pixel. The Optimize Cameras tool could then refine the first alignment by optimizing the camera positioning on the remaining points. The results of the first cloud-to-cloud alignment operated on separated chunks of the same site is a set of 4 different chunks, one for each study area. Dense cloud parameters were set on Medium quality and Aggressive Filtering to remove small outliers. The number of points generated in the dense cloud process (Table 3) depends on the site area and image overlap, representing the paramount quality of metrical data produced during the investigation. Even high- quality meshes could not fit this amount of information or, at least, no computer could manage operations with similar numbers of polygons. For this reason, joint measurement was performed on the dense clouds instead of meshes.

4.2.2. Generation of 3D Models

2.5D and 3D models were generated creating a polygon structure around the points provided by the dense clouds. To achieve a higher resolution model, 2.5D reconstruction was driven instead of a 3D one, reducing computation time and providing a high-quality mesh of vast areas. The difference between 2.5D and 3D meshes, both classified as 3D models, lies in the single vertical coordinate for a point in a 2.5D mesh. Although perfect vertical structures lose their detail and undercuts, the horizontal structures gain higher quality in a comparison of the same face number 3D model. Using the dense cloud for measurement and the 3D only for DEM and orthomosaic generation, no data are lost in the procedure. High Field value was set instead of Arbitrary to generate a 2.5D model, the base of the construction was the dense cloud and the face count was manually set according to the area dimension and complexity in a range from 20 to 30 million polygons. After 3D modeling of the entire workspace, DEMs were generated and exported according to the WGS84 reference system.

Before texture processing, another picture selection was performed to improve the quality of the results, dealing with high redundancy of data. If for alignment and dense clouds generation every picture with 0.70 or higher quality value was employed to avoid blurry areas in the texture a refinement selection was performed. Every picture with a quality evaluation lower than 0.80 was deselect in order to avoid their contribution in the texturing process. Mosaic texture mapping performs well with high-quality input files avoiding the low-quality loss related to the ‘Average’ method.

With the same picture selection, the orthomosaics were assembled with geographic projection and the same blending mode of the texture. “Enable hole filling” was selected even if a 2.5D model usually does not have holes except on the edges of the model. The pixel size and dimension value were automatically set at the higher value permitted by the original picture resolution.

4.3. Gravity-Induced Joints and Coastal Megaclast Deposit Investigation

Analyses of orthomosaics and 3D models permitted the identification and property investigation of the main discontinuities originated or enlarged by rock spreads. In particular, we analyzed the persistence of joints and variations of apertures of each joint every 5–10 m, not exceeding distances of 10 m. Persistence is the term used to describe the extent or size of a joint within a plane. It has a major influence on the fragmentation characteristics during detachment and fall of external portion of limestone plateau during lateral spreading and block sliding processes. Table 4 lists the categories to classify joint persistence used in this paper according to the ISRM Commission [76].

3D models were also used for characterizing the changes in joint apertures every 5–10 m along their length. Joint aperture is the normal distance separating the adjacent walls of open discontinuity. A set of terms to define aperture values was suggested by ISRM [76]. Table 5 lists the aperture terms and associated width.

Finally, we used orthomosaics and 3D models for the detection of blocks and their size categorizations using the three classes proposed by Bruno and Ruban [55]. It is important to emphasize that the measurement of the major axis of several block was realized using 3D models, which permit a better understanding of the real geometries of slope-failure accumulations. We used 3D models that are fundamental for the GIS operator in the recognition of the blocks by visualizing in a stereoscopic environment avoiding misinterpretation due to 2D visualization, as occur using GE images. The results achieved using images provided by GE in terms of geometry and size are listed in the next paragraph, whereas UAV-DP outputs are successively discussed.

5. Results of GE Analysis

The GE analysis carried out on block slides’ megaclast deposits enabled the identification and categorization of a total population of 5059 megaclasts, which were digitized and stored in a QGIS environment. The GE imagery used in the elaborations was retrieved from satellite Landsat 8, acquired on August 2nd, 2019 (the more recent data available at the time of this research), with a resolution of 15 m. Figure 8 shows their spatial distribution in the four study sites. Table 6 lists the category ranks of each study area.

Although the block size category dominates other size classes in all four localities, megablocks are abundant in Il-Qarraba peninsula and Ras-Il-Wahx headland where there are also four superblocks. A superblock located along the northern side reaches almost 300 m.

6. Results of UAV-DP Analysis

The main outcomes of DP technique were orthomosaics and 3D models. These products permitted the identification, mapping, and categorization of gravity-induced joints located along external portions of UCL plateaus and blocks making up large coastal accumulations situated between the coastline and cliff escarpments. In addition, UAV-DP allowed the production of inventories of both gravity-derived joints and megaclast deposits, described and commented as follows.

6.1. Joint Inventory and Property Characterization

Joints originating from rock spreads are clearly visible in DP outputs computed using Agisoft Photoscan^TM, related to the fact that these discontinuities are often ‘cavernous’ joints and their lengths can exceed 100 m.

We categorized the tension cracks as (1) single joint, and (2) complex network of discontinuities. Single joints are generally internal to the plateaus and widen as a result of lateral spreading. A complex network of discontinuities is intended as a dense net of open and vertical joints located at the edge of the plateaus. These discontinuities are laterally unconfined. Single joints indicate an initial or intermediate stage of lateral spreading, whereas a dense network of cracks corresponds to an advanced stage of lateral spreading that often evolve into block sliding, as evidenced by the high number of blocks located downslope and covering BC. Figure 9 shows the location of single joints and complex networks of discontinuities inventoried using orthomosaics and 3D models obtained by UAV-DP. The base maps used are orhophotos produced by the UAV-DP technique.

We identified seven main single persistent discontinuities. Except for Joint #7, they are all characterized by ‘very high persistence’ (>20 m). Table 7 lists their persistence values and maximum apertures.

The inner sector of Anchor Bay is the area most populated by single discontinuities (Figure 9B). Their occurrence proves an intermediate stage of lateral spreading confirmed by active deformations recorded by GNSS monitoring [59,60]. We used 3D models to determine variations in aperture along the length of joints every 5–10 m. Figure 10 lists the aperture values of seven joints associated with the purple points listed in Figure 9.

Aperture values are important indicators of the lateral spreading stages. Low aperture values correspond to a recent formation, whereas high values correspond to an intermediate stage of lateral spreading. According to ISRM ranks [76], Joint #3 already described in paragraph 3.3 is a ‘cavernous discontinuity’ characterized by extremely high values of aperture exceeding 3 m. This ‘very high persistence‘ joint is accompanied by Joint #4 and Joint #5, which are located easterly and are characterized by lower values of apertures with respect to Joint #3.

The first well-developed network of vertical discontinuities is located in the westernmost part of Marfa Ridge limestone plateau and is made up of 11 discontinuities (Figure 11A). Figure 11 also shows a detailed map of a second complex network of discontinuities (17 joints) detected at Il-Qarraba (Figure 11C) and an oblique photo of joints taken by a UAV flight (Figure 11B).

The charts in Figure 12 correlate the persistence (x-axis) and maximum values of apertures (y-axis) of two complex networks of discontinuities.

The joints of these two complex networks are often cavernous, vertical, and tens of meters deep, showing an advanced phase of rock spreading. The total population of Marfa Ridge joints is included in the ‘cavernous’ class or ‘extremely wide’ class. A similar condition occurs at Il-Qarraba caprock, where only two joints are tight. These advanced phases of displacement are confirmed by the high population of blocks downslope these networks, and the field observations.

The megablocks located in front of limestone cliffs are rotated and lowered, showing a modification of the landslide mechanism from lateral spreading to block sliding. The rotation and lowering movements of wide megablocks are clearly visible in the front of Il-Qarraba and Marfa Ridge. These motions are often contemporary and are consequences of a vertical motion related to the occurrence of a slip surface. These slip surfaces are steep along the external edge of the UCL plateaus whereas their angles become almost planar in the areas adjacent to the coastline [39].

6.2. Coastal Megaclast Deposit Inventory and Categorization

Orthomosaics and 3D models produced by DP permitted the counting and categorization of 8943 megaclasts. Figure 13 shows the spatial distribution of these megaclasts categorized according to the methodology developed by Bruno and Ruban [55]. Blue and yellow points indicate block and megablock populations, respectively. The wide mustard-yellow areas located in Figure 13C represent four superblocks. Table 8 lists the category ranks of each study site.

Although the block size category dominates other two classes in all the study sites, the megablock size population is high mainly in Il-Qarraba peninsula and Ras-Il-Wahx headland (Table 8), where there are also four superblocks.

The low population of megaclasts in Anchor Bay is related to the small study area and to the intermittent occurrence of blocks on the inner side of the cove, which is more protected by active sea erosion than the external area [60].

7. Discussion

This research permitted the production of a detailed inventory of gravity-induced joints and megaclast deposits of the northwestern coast of Malta, showing the benefit of using micro drones in the study of large slow-moving landslides. In particular, lateral spreads and block slides are abundant in this sector of the island and were recently investigated by performing geomorphological mapping assisted by visual interpretation of aerial and satellite images, as reported in the literature [14,39,40,49]. The latter provided detailed data regarding their spatial distribution, including the characteristics, size, and their main geomorphological features.

The drawback of conventional geomorphological field mapping in block sliding areas is represented by the limited accessibility of sites to be surveyed, due to hummocky terrains often dangerous for field operators. For this reason, visual interpretation of photographs taken at low-altitude can facilitate geomorphological investigations; micro drones represent a valid and inexpensive instrument to support them. Micro drones speed up the surveys and reduce the cost of operation activities in the field. In addition, their images show higher resolution compared to the GE images. Micro drones are excellent solutions because of their portability, lightness, not bulky and versatility in impervious landslide areas, such as those of northwestern Malta. In addition, commercial drones with weights around 0.3 and 1 kg remain the most inexpensive type of UAVs, as reported by Hassanalian and Abdelkefi [70].

The Maltese UAV flight regulations and the rather dry climate of Malta facilitate the use of UAVs and avoid the drawbacks listed by Casagli et al. [77], such as restrictive flight protocols and the negative effect of vegetation on the point cloud. Flight regulations of Malta are very clear whereas vegetation cover is rare and limited to low-growing scrub [78].

Concerning gravity-induced joints related to lateral spreading occurring along the northwestern coast of Malta, a suitable agreement was found between existing datasets and UAV-DP outputs in terms of the identification of joints and their persistence (Figure 14). This is due to the clear difference in color between limestone pavement (light gray) and each joint. The latter is often deep and is witnessed by a black color, which facilitates the identification and mapping also using GE images.

Conversely, the use of the UAV-DP technique shows great advantages for joint aperture determination and distribution outline (Figure 15). Differences in aperture values can exceed 1 m, and are often related to the few field aperture measurements, which are a highly time-consuming task along ‘very high persistence’ joints. It should be noted that Figure 15 shows at comparison of 17 joints (out of 35) at Marfa Ridge and Il-Qarraba due to the fact that some joints were not identified during field surveys or aperture measurements in the field were not carried out for safety reasons (joints located on unstable rock masses).

In general, the joint aperture is very sensitive to rapid changes in loading conditions or deviations of stress scenarios [79] and thus its value can vary sensibly along discontinuities exceeding 100 m in length, as in Anchor Bay (Joint #3) or Marfa Ridge (Joint #1). The UAV-DP technique and its products are capable of measuring aperture changes of each tension crack remotely, avoiding a time-consuming part of field activities related to joint investigations. A total of 132 aperture values were measured in 3D models and stored in a QGIS environment.

The large amount of accurate joint property data collected for this research may constitute a survey “0”, as it may be compared to UAV-DP outputs obtained by the same flight plans stored in Litchi platform. The temporal analysis of the variation of joint properties represents an important step for a better understanding of lateral spreading evolution.

With reference to coastal megaclast deposits, the GE analysis carried out for this research highlighted thousands of megaclasts widespread at the toe of the limestone cliffs (Table 6). Their spatial distribution shows areas with a low occurrence of megaclasts, related to the sector of cliffs with a low density of fractures (Figure 8). The ratio between blocks and megablocks varies between 9.68 at Ras Il-Wax and 21.29 at Marfa Ridge. This ratio increases using the DP and ranges between 9.3 at Ras Il-Wax and 38.28 at Marfa Ridge. These differences are mainly due to the higher resolution of DP-derived orthomosaics.

Significant differences between the results of UAV-DP outputs and GE analysis are evident in the counting of megaclast populations in the four localities, as reported in Figure 16.

The differences of megaclasts observed at each location are remarkable, mainly at Marfa Ridge and Il-Qarraba peninsula. The visual interpretation of GE images did not allow a detailed characterization of megaclasts surrounding the limestone plateaus.

As reported in a recent study of four localities populated by coastal megaclast deposits [80], GE image resolution does not permit a clear distinction between megaclasts, and their identification is not always easy, especially when the slope is irregular. Conversely, megaclasts displaced by block slides are easily distinguished using 3D models or orhomosaics obtained by low-altitude drone flights.

The recognition of megaclasts using GE can be significantly reduced by the shadow effects at the cliff toe. A clear example is the insufficient visibility related to the shadow effect of the vertical cliffs that limit the plateau. The steep cliffs generate wide shadowy areas at the cliff toe that mask most of the 1–10 m blocks, preventing their recognition. Unfortunately, image post-processing does not help. Figure 17A displays an example of a view from a HR image provided by GE with a 14 m wide shadow area at Marfa Ridge plateau. The same area is well lit in Figure 17B, having chosen the most suitable moments for the drone survey.

8. Conclusions

Research on lateral spreads and block slides normally requires long-term investigations due to their ‘extremely slow’ or ‘very slow’ geomorphological evolution. Their analysis needs extensive field surveys, which are time-consuming and expensive tasks. For these reasons, the recent development of commercial drones weighting less than 1 kg can be a valid support for field surveys. Micro drones can be used to rapidly map large areas and provide VHR images. UAV-DP can integrate field surveys with reference to the detection and property characterization of ‘very high persistence’ joints, and significantly improve the quality of block sliding investigations, in particular with respect to satellite image analysis using GE.

Gravity-induced joints are typically rough and extensive, located on the edge of plateaus and are crucial for the recognition of the internal perimeter of areas affected by rock spreads. We identified 35 main joints in four study sites, most of which grouped into two complex networks of discontinuities. These joints range between 5 m and 158 m of Joint #1 located on Marfa Ridge. Most of these joints are characterized by high aperture values, which can exceed 3 m, such as for Joint #3 located in the central sector of Anchor Bay. UAV-derived investigations of fracture aperture are very important, especially when discontinuities limit blocks that are inaccessible or highly unstable and therefore dangerous [81]. The application of the UAV-DP technique permitted to accurately estimate aperture values of 18 joints neither identified on the field nor surveyed for safety reasons, and to improve the number of aperture measurements along ‘very high persistence‘ joints, with a great increase of aperture data of 600% compared to field surveys (132 observations vs. 22 field measurements). This improvement is due to the possibility of investigating and measuring the aperture values of joints using 3D models, which avoids field visits along joints and rock mechanic measurements that are time-consuming and unsafe.

Regarding the extensive megaclast deposits that populate the northwestern coast of Malta island, UAV-DP products permitted the identification and categorization of 8576 blocks (1–10 m), 363 megablocks (10–100 m) and 4 superblocks (over 100 m and all surrounding Ras-Il Wahx). The comparison with GE data shows a substantial difference with respect to UAV-DP investigations. The total population of megaclasts in the four localities varied from 5059 individuals identified with GE to a new DP-derived population of 8943 individuals, with an increase in the total population of 76%. The principal differences were recorded at Il-Qarraba and Marfa Ridge, where an advanced phase of lateral spreading and block sliding can be inferred. This evolution assumption of slope-failure processes can be confirmed by extremely low rates of surface displacements recorded at Il-Qarraba by repeated GNSS surveys.

The movement of the blocks on clayey slopes has increasingly fragmented the megaclasts on their way towards the sea, decreasing the volume of many megaclasts, and thus making it difficult for GIS operators to recognize them. In addition, the low visibility of clayey slopes in GE images, caused by the steep plateau cliffs, can explain the large population differences between GE and UAV-DP outputs recorded at Marfa Ridge and Il-Qarraba peninsula.

Quantitative information about tension cracks and megaclast deposits collected from the DP technique is an important step for a better understanding of large landslides such as rock spreads and block slides. The inexpensive technique (approximately 2k USD for drones plus 3k USD for software) used in this study shows the potential for being extended to other rocky coastal sites affected by these types of slow-moving landslides.

Further developments of this research will include the collection and analysis of the new aerial photos taken within the same flight plans to identify the occurrence of new or enlarging joints, and the measurement of heights of steep cliffs of UCL plateaus using 3D models. These values will be more accurate than cliff heights obtained by Maltese base maps and might be used for slope stability analysis.

Author Contributions

All authors have contributed to this study. S.D. conceived the paper; S.D., and V.M. performed the UAV flights; V.M. computed DP; Original idea: S.D., S.F., and M.M.; Conceptualization: all authors; Writing—Original draft: all authors; Supervision: M.S.; Writing—Review & editing: S.D., and M.S. All authors have read and agreed to the published version of the manuscript.

Funding

This project is part of the Project “Coastal risk assessment and mapping” funded by the EUR-OPA Major Hazards Agreement of the Council of Europe (2020–2021). Grant Number: GA/2020/06 n_654503 (Resp. of University of Modena and Reggio Emilia Unit: Mauro Soldati).

Acknowledgments

We thank Linley Hastewell (University of Portsmouth) for his comments and suggestions during the preparation of the manuscript. We are also grateful to the three reviewers for their useful comments.

Conflicts of Interest

The authors declare no conflict of interest. The funders and DJI^TM (Nanchan District, Shenzhen, China) had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the paper, or in the decision to publish the results.

Abbreviations

The following abbreviations are used in this paper:

GE	Google Earth
UAV	Uncrewed Aerial Vehicle
DP	Digital Photogrammetry
DGSD	Deep-seated Gravitational Slope Deformation
RPAS	Remotely Piloted Aircraft System
HR	High Resolution
VHR	Very High Resolution
DEMs	Digital Elevation Models
LiDAR	Light Detection And Ranging
ALS	Airborne Laser Scanning
TLS	Terrestrial Laser Scanning
LIM	Landslide Inventory Map
UCL	Upper Coralline Limestone Formation
BC	Blue Clays
GL	Globigerina Limestone Formation
LCL	Lower Coralline Limestone Formation
GNSS	Global Navigation Satellite System
A-DinSAR	Advanced SAR interferometric analysis
FS	Fast Static
UAS	Uncrewed Aircraft System
VL/MC	Very Light MultiCopter
MALE	Medium-Altitude Long-Endurance
HALE	High-Altitude Long-Endurance
VTOL	Vertical TakeOff and Landing

Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Figures and Tables

Figure 1. (A) Location of Malta; (B) Map of the Maltese islands; and (C) Relief map of the northwestern sector of Malta.

Figure 2. Landslide distribution along the northwestern coast of the island of Malta.

Figure 3. Landslide evolution along the northwestern coast of Malta. Schematic and simplified examples of (A) early stage of rock spreading; (B) lateral spreading evolving into block sliding; and (C) advanced stage of landsliding (cf. [52]).

Figure 4. Views of the study sites: (A) Marfa Ridge; (B) Anchor Bay; (C) Ras Il-Wahx headland; and (D) Il-Qarraba peninsula.

Figure 4. Views of the study sites: (A) Marfa Ridge; (B) Anchor Bay; (C) Ras Il-Wahx headland; and (D) Il-Qarraba peninsula.

Figure 5. Map showing the main landforms of Il-Prajjet (background GE image). White arrows point to the open joint, which partially isolates and displaces a 43,500 m3 rock slice. The light grey areas highlight downslope megablocks.

Figure 6. Map showing the location of Il-Qarraba caprock and its rock masses displaced by rock spreading and block sliding. Tilted megablocks transported by block slides are clearly visible (background GE image).

Figure 7. Workflow of GE imagery analysis and UAV-DP procedures.

Figure 8. The spatial distribution of categorized megaclasts of: (A) Marfa Ridge; (B) Anchor Bay; (C) Ras Il-Wahx headland; and (D) Il-Qarraba peninsula.

Figure 9. Joint maps obtained by 3D-model analyses of: (A) Marfa Ridge; (B) Anchor Bay; (C) Ras Il-Wahx; and (D) Il-Qarraba. Numbers inside yellow squares indicate single points. The purple numbers located along the joints approximately indicate the locations where aperture observations were performed using 3D models.

Figure 10. DP analysis of aperture variations of: (A) Joint#1 and Joint#2; (B) Joints located in the central and inner part of Anchor Bay; and (C) Joints located along the edge of Ras Il-Wahx plateau.

Figure 11. (A) Oblique view of complex network of discontinuities located in the western part of Marfa Ridge; (B) Spot on the western part of Il-Qarraba caprock where Rock mass #II, Rock mass #III, and Rock mass #IV are partially isolated by persistent joints; and (C) UAV-DP derived map of complex network of discontinuities.

Figure 12. Analysis and ISRM ranking of persistence and aperture of joints belonging to complex network of discontinuities of (A) Marfa Ridge; and (B) Il-Qarraba caprock.

Figure 12. Analysis and ISRM ranking of persistence and aperture of joints belonging to complex network of discontinuities of (A) Marfa Ridge; and (B) Il-Qarraba caprock.

Figure 13. Megaclast inventory maps of (A) Marfa Ridge; (B) Anchor Bay; (C) Ras Il-Wahx headland; (D) Il-Qarraba peninsula; (E) percentages of each classes with respect to the total population.

Figure 14. Averaged length per joint using GE-derived existing dataset and UAV-DP outputs. The persistence terms follow ISRM standards [76].

Figure 15. Joint apertures determined using field measurements from existing dataset and DP-derived 3D models in the four localities.

Figure 16. Megaclast populations according to the methodology proposed by Bruno and Ruban [55] using GE analysis and UAV-DP outputs.

Figure 17. Identification and size categorization of Marfa Ridge megaclast deposit using (A) GE image where a 14 m wide shadow area complicates the recognition of megaclasts at the foot of the vertical escarpment that bounds the plateau; and (B) UAV-DP derived orthomosaic. The absence of an extensive shadow assists the GIS operator to identify in the area outlined by the dashed rectangle a dozen of blocks not recognized using GE image.

Table 1

Main characteristics of the DJI Mavic^TM.

Weight[kg]	Size[mm]	Max Speed[km/h]	Flight Time[min]	Ground Resolution(100 m–48 mm)	Max Mapped Surface(1.9 m/pix–80% Overlap)
0.905	322 × 242 × 84	72	~30	0.019 m/pix	0.1 km2

Table 2

Classification of UAVs following the weight methodology of Brooke-Holland [68].

Type	Weight Range [kg]
Nano drones	0–0.2
Micro drones	0.2–2
Mini drones	2–20
Small drones	20–150
Tactical drones	150–1500
MALE, HALE, strike drones	>1500

Table 3

Information on flights performed in the study sites.

Site	Flight Time [min]	Picture [#]	AreaCoverage[m2]	Covered Flight Distance [km]	Dense Cloud Model Point Count	3D ModelFace Count
Marfa Ridge	37	838	698,769	12.2	90,370,838	18,367,407
Anchor Bay	13	332	246,980	3.6	64,553,718	12,910,687
Ras Il-Wahx	69	764	671,447	18.6	189,249,228	11,667,442
Il-Qarraba	36	886	267,656	9.7	142,510,544	9,464,093

Table 4

Classification of joint persistence [76].

Type	Joint Length [m]
Very low persistence	<1
Low persistence	1–3
Medium persistence	3–10
High persistence	10–20
Very high persistence	>20

Table 5

Classification of joint aperture [76].

Description	Width of Aperture
Very tight	<0.1 mm
Tight	0.1–0.25 mm
Partly open	0.25–0.5 mm
Open	0.5–2.5 mm
Moderately wide	2.5–10 mm
Wide	Over 10 mm
Very wide	1–10 cm
Extremely wide	10 cm–1 m
Cavernous	>1 m

Table 6

Megaclast size analysis within the study sites using GE images.

Location	Block	Megablock	Superblock
Marfa Ridge	1448	68	-
Anchor Bay	479	23	-
Ras Il-Wahx	968	100	4
Il-Qarraba	1861	108	-

Table 7

Main joints recognized by 3D models analysis.

Joint [#]	Location	Orientation	Persistence [m]	Max Aperture [m]
1	Marfa Ridge	NW-SE	158	0.75
2	Anchor Bay	NW-SE	86	1.14
3	Anchor Bay	W-E	109	3.25
4	Anchor Bay	W-E	53	0.49
5	Anchor Bay	W-E	34	0.98
6	Ras Il-Wahx	NW-SE	57	2.03
7	Ras Il-Wahx	NW-SE	16	0.64

Table 8

Megaclast size analysis within the study sites using drone DP outputs.

Location	Block	Megablock	Superblock
Marfa Ridge	3369	88	-
Anchor Bay	760	42	-
Ras Il-Wahx	1008	108	4
Il-Qarraba	3439	125	-