**Sea Level Rise Scenario for 2100 A.D. in the Heritage Site of Pyrgi (Santa Severa, Italy)**

**Marco Anzidei 1,\*, Fawzi Doumaz 1, Antonio Vecchio 2,3, Enrico Serpelloni 1, Luca Pizzimenti 1, Riccardo Civico 1, Michele Greco 4, Giovanni Martino <sup>4</sup> and Flavio Enei <sup>5</sup>**


Received: 29 November 2019; Accepted: 15 January 2020; Published: 21 January 2020

**Abstract:** Sea level rise is one of the main risk factors for the preservation of cultural heritage sites located along the coasts of the Mediterranean basin. Coastal retreat, erosion, and storm surges are posing serious threats to archaeological and historical structures built along the coastal zones of this region. In order to assess the coastal changes by the end of 2100 under the expected sea level rise of about 1 m, we need a detailed determination of the current coastline position based on high resolution Digital Surface Models (DSM). This paper focuses on the use of very high-resolution Unmanned Aerial Vehicles (UAV) imagery for the generation of ultra-high-resolution mapping of the coastal archaeological area of Pyrgi, Italy, which is located near Rome. The processing of the UAV imagery resulted in the generation of a DSM and an orthophoto with an accuracy of 1.94 cm/pixel. The integration of topographic data with two sea level rise projections in the Intergovernmental Panel on Climate Change (IPCC) AR5 2.6 and 8.5 climatic scenarios for this area of the Mediterranean are used to map sea level rise scenarios for 2050 and 2100. The effects of the Vertical Land Motion (VLM) as estimated from two nearby continuous Global Navigation Satellite System (GNSS) stations located as close as possible to the coastline are included in the analysis. Relative sea level rise projections provide values at 0.30 ± 0.15 cm by 2050 and 0.56 ± 0.22 cm by 2100 for the IPCC AR5 8.5 scenarios and at 0.13 ± 0.05 cm by 2050 and 0.17 ± 0.22 cm by 2100, for the IPCC Fifth Assessment Report (AR5) 2.6 scenario. These values of rise correspond to a potential beach loss between 12.6% and 23.5% in 2100 for Representative Concentration Pathway (RCP) 2.6 and 8.5 scenarios, respectively, while, during the highest tides, the beach will be provisionally reduced by up to 46.4%. In higher sea level positions and storm surge conditions, the expected maximum wave run up for return time of 1 and 100 years is at 3.37 m and 5.76 m, respectively, which is capable to exceed the local dune system. With these sea level rise scenarios, Pyrgi with its nearby Etruscan temples and the medieval castle of Santa Severa will be exposed to high risk of marine flooding, especially during storm surges. Our scenarios show that suitable adaptation and protection strategies are required.

**Keywords:** sea level rise; coastlines; 2100; storm surges; heritage sites; Pyrgi; Mediterranean; UAV; DSM

#### **1. Introduction**

Observational and instrumental data collected worldwide since the last two-three centuries show that the global sea level is continuously rising with an accelerated trend in recent years, which coincides with the rise in global temperatures [1]. Global mean sea level is expected to rise by about 75 to 200 cm by 2100 in the worst scenarios [2–5], i.e., the most serious effects of climate change that might occur in future decades. These values will be even larger in subsiding coasts of the Mediterranean, entailing widespread environmental changes, coastal retreat, marine flooding, and loss of land, which will be disadvantages for human activities. The sea level rise will amplify the impacts exerted by a multitude of hazards (i.e., storm surges, flooding, coastal erosion, and tsunamis) on the infrastructure and building integrity, people safety, economic assets, and cultural heritage.

Therefore, it is important to mitigate these risks by providing multi-temporal scenarios of expected inland extension of marine flooding as a consequence of the sea level rise for a cognizant coastal management [6–8].

This is particularly true for the Mediterranean region where ancient civilizations were born and developed along its coasts [9,10]. A large number of heritage sites are located at the waterfront or very close to the sea level and are exposed to marine flooding under the effects of ongoing climate change. A large part of these sites, which are dated back to the Greek, Roman, and Medieval ages, are exposed at increasing risks to coastal hazards that are related to a sea-level rise [10].

Aerial photogrammetric surveys performed by small Unmanned Aerial Vehicles (UAVs) can provide accurate topography at low costs and, in short time, for small areas with respect to conventional aerial surveys [11]. Results, when analyzed in combination with sea level rise projections and vertical land movements (VLM), can support the realization of the expected sea-level rise and storm surge scenarios for future decades [12].

In this study, we show an effective application for the coastal archaeological area of Pyrgi, Italy, which is near Rome (Figure 1). High resolution maps of expected flooded areas and coastal positions for 2100, even in storm surge conditions, are reported in this paper.

**Figure 1.** The investigated area of Pyrgi with the location of the GNSS stations of MAR8 and TOLF and the nearby tide gauge station of Civitavecchia.

Our approach is to apply a multidisciplinary methodology previously tuned in the savemedcoasts project (www.savemedcoasts.eu), which includes topography, geodesy, sea level data, and climatic projections to estimate realistic sea level rise scenarios for targeted coastal areas. Our approach provides a useful analytical tool to identify the best adaptation and defense strategies against the sea level rise impact, and to protect heritage sites. The results help decision-makers in the selection of the best practical actions aimed at preserving the archaeological and historical sites located in coastal areas that are subjected to sea level-related risks. The proposed methodology can be exported in other areas of the Mediterranean region and beyond its borders.

#### **2. Geo-Archaeological Setting**

The heritage site of Pyrgi is located along the coasts of Northern Latium, between the villages of Santa Severa and Cerveteri, which is about 50 km north of Rome (Italy) (Figure 1). The area includes the Castle of Santa Severa that, with Pyrgi, is one of the most important heritage sites of the Tyrrhenian coast. The area has been settled since the V-IV millennium B.C. [13] and continuously developed during the Neolithic Age, during the Bronze Age (II millennium B.C.), and during the Iron Age (IX-VIII century B.C.), thanks to its good environmental conditions. In the Etruscan phase (VII-IV century B.C.), Pyrgi was the port of the ancient Etruscan city of Caere and played an important role in the maritime commerce being frequented by Greeks and Phoenicians ships. The area includes a sanctuary and the temples of Eileithyia-Leukothea and Apollo, Cavatha, Suri, and the Etruscan Uni analogue to the Phoenician Astarte [13].

After the Romanization of this area (III century B.C.), Pyrgi became a maritime colony and a tall fortress surrounded by a polygonal wall was built on part of the Etruscan settlement.

During the Roman imperial age, the city of Pyrgi continued to be frequented until the 5th-6th century A.D. when a Byzantine castrum was built on its remains with the early Christian church of Santa Severa inside.

Later, the medieval and renaissance village became a large farm located in a strategic position between the main harbors of Rome and Civitavecchia [13–16]. In terms of its geological setting, the coast of Pyrgi belongs to the roman co-magmatic province [17] that underwent major volcanic activity during the Plio-Pleistocene era (Figure 2). The surface geology is characterized by a bedrock belonging to the allochthonous Outer Ligurids [18], represented by the Tolfa Formation, spanning from the Late Cretaceous to Palaeogene [19], defined as the Pietraforte formation and "comprehensive succession" [20]. The latter consists of a series of marl limestone and grey marl beds that outcrop between Rome and Civitavecchia underlying the Pietraforte unit. Biogenic sandstones of the Early-Middle Pleistocene Panchina Formation partly overlies the Pietraforte [21].

**Figure 2.** Simplified geological map of the study area (from http://dati.lazio.it/catalog/it/dataset/cartageologica-informatizzata-regione-lazio-25000). Legend: (1) anthropic debris, (2) coastal and marsh sands and recent dunes, (3) calcareous marls and clays, (4) flysch, (5) clays with chalks, (6) landslide de posits, (7) alluvial deposits, (8) lavas, (9) travertines, (10) sandy deposits, and (11) sandy deposits of marine facies.

The Holocene deposits are represented by travertines, slope debris, alluvial, weathering deposits, gravels, and sandy beaches [22]. The neotectonic of this sector of the Tyrrhenian coast of Italy is marked by the elevation of the MIS 5.5 marine terraces that show stability and a weak uplift of the inland sector for the last 124 years [23,24], which is related to magmatic injections under the Vulsini and Sabatini volcanic complexes [24]. Reference [25] underlines that the long-term uplift may not be an appropriate description for all the past two millennia since some weak subsidence may have occurred at Pyrgi and along the nearby coasts.

#### **3. Methods**

We applied a multidisciplinary approach using coastal topography, geodesy, and climatic-driven estimates of the sea-level rise to provide maps of flooding scenarios for the year 2100 A.D. for the coast of Pyrgi. Our study consists of three main steps: (1) the realization of UAV surveys to obtain an ultra-high-resolution orthophoto and a DSM model of the coastal area to map the current and the projected coastline positions including sea level data in the analysis, (2) the estimation of the current vertical land movements from the analysis of geodetic data from the nearest GPS stations, and, (3) by combining these data with the regional IPCC-AR5 projections (RCP-2.6 and RCP-8.5 scenarios), the calculation of the upper bounds of the expected sea levels for the targeted epochs of 2050 and 2100 A.D. and the corresponding expected inland extent of the marine flooding and shoreline positions were calculated. Lastly, storm surge scenarios were implemented for return times of 1 and 100 years, for sea level rise conditions.

#### **4. Digital Terrain Model Reconstruction**

To realize the ultra-high-resolution Digital Surface Model (DSM), an aerial photogrammetric survey was performed using a radio-controlled multi-rotor Da Jiang Innovation (DJI) Phantom 4pro UAV system, equipped with a high resolution lightweight digital camera, to capture a set of aerial images of the investigated area (Table 1).


**Table 1.** Survey features.

The UAV was controlled by an autopilot system using waypoints previously planned by PIX4D® capture IOs App as a Ground Control Station system. To optimize the photogrammetric spatial resolution and coverage of the surveyed area, a constant altitude of 70 m was maintained during the flight and 306 partly overlapping (70% of longitudinal and 70% lateral overlapping between subsequent photos) aerial digital photos were acquired during three successful flights of 13, 11 and 15 min of duration each. To scale the aerial images, we used a dual-frequency geodetic Global Positioning System Real Time Kinematic (GPS/RTK) receiver STONEX S900A® to measure the coordinates of a set of reference Ground Control Points (GCPs) falling in the investigated areas. GCPs positions were estimated in real time by the RTK technique with about 1 to 2 cm of accuracy, with respect to the reference GPS station TOLF. The latter is part of the GNSS network of the Istituto Nazionale di Geofisica e Vulcanologia (INGV) [26] and also pertains to the Leica SmartNet ItalPoS network (https://hxgnsmartnet.com/en-gb/) for real-time positioning services. We used the Agisoft PhotoScan® software package (http://www.agisoft.com) based on the Structure-from-Motion photogrammetry technique [27] to process the acquired georeferenced images. The analysis included: (1) camera alignment with image position and orientation, (2) generation of a dense points cloud, and (3) generation of an orthophoto covering a land surface of 0.226 km2 with a ground resolution of 1.94 cm/pixel as well as the DSM creation in the WGS84-UTM32 coordinate system.

The obtained ortho-rectified images (orthophotos) and digital elevation models were also managed by Global Mapper® and ESRI ArcGis® software to create cross sections, slope maps, surfaces, and coastline positions as well as calculate the dimension of the potential flooded areas. The extracted Digital Surface Model has a resolution of 15.5 cm/pixel and a point density of 41.5 points/m2

As targets for GCPs measured by GPS/RTK, we used a set of i) natural markers belonging to fixed structures (e.g., the center of manholes, wall and sidewalk corners, small structures, or large stones), and ii) mobile targets for the time of the flight, such as thin metal crosses of 60 × 60 cm of size, deployed in the investigated area. All these GCPs were chosen as areas recognizable on the images during data analysis.

In total, 22 Ground Control Points (GCPs) and Check Points (CPs) were used to geo-reference the orthophoto with the ground control point toolset of AGISOFT photo scan (Figure 3). The 3D coordinates of these points have been estimated with a mean RMS of 0.6 and 0.9 cm for the horizontal and vertical components, respectively, and were used to evaluate the vertical accuracy of our final DSM. Figure 3 shows the orthophoto while Figure 4 shows the DSM.

**Figure 3.** The orthophoto with the Ground Control Points (red dots) position used during the UAV surveys.

**Figure 4.** The Digital Surface Model (DSM) of the coastal sector of Pyrgi, from the analysis of the aerial photos.

#### **5. Tidal Correction and Coastline Position**

The coast of Pyrgi, similar to most of the coasts of the Mediterranean Sea, is characterized by a micro-tidal environment and tides are generally in the range of ±30 cm. Only in the Gulf of Gabes (Tunisia) and the North Adriatic Sea (Italy), tides may reach amplitudes up to about2m[27–30]. We used the tidal data collected by the Italian tidal network managed by the Italian Institute for Environmental Protection and Research (ISPRA) (data are freely available at www.mareografico.it) at the sea level station of Civitavecchia (located at LAT 42◦05 38.25", LON 11◦47 22.73 ), which is placed near Pyrgi (Figure 1), to estimate the tide level (TL) and the local mean sea level (LMSL) at the time of the UAV surveys (26 April, 2019 at 07:30 or 05:30 UTC). We considered the complete time series to account for a long-term linear trend, representative of the mean sea level during UAV surveys, with respect to which the TL is defined. This tide gauge station shows a valid recording period of about eight years (2011–2019) with a sea level trend of 0.25 ± 0.1 mm/a, which is calculated from a linear fit on the monthly data (Figure 5a). From the data analysis, a mean tide amplitude of ~35 cm and a maximum tidal range up to ~60 cm has been estimated (Figure 5b). To define a reference level for elevation data, the mean sea level was computed propagating the linear trend from the time of surveys, assuming the mean sea level as a reference value for the year 2018. A mean sea level of 4.8 ± 11.8 cm above the topographic benchmark was estimated from the hourly tidal data. Since the UAV surveys have been performed during a high tide of 4 cm, as shown by the tidal recordings (www.mareografico.it), then the reference sea level at the time of surveys corresponds to a position of only 0.8 cm above the mean sea level for 2018.

This value is negligible compared to the DSM accuracy, so we did not consider it for the analysis of sea level rise scenarios. Lastly, we used the local mean sea level calculated at the tide station of Civitavecchia as elevation data, given the small distance between this location and Pyrgi. The obtained value was used to define the position of the coastline during the surveys and the one expected from the sea level rise scenarios to 2100 A.D.

**Figure 5.** Sea level data analysis for the tide gauge of Civitavecchia, which is located a few km north of Pyrgi (see Figure 1 for location). (**a**) Monthly data of sea level recordings collected in the time span 2011–2019 (about nine years). The red line is the linear fit of the sea level trend at 0.25 ± 0.1 mm/yr. (**b**) Statistical diagram of sea level heights (cm) versus time (hours/year). The values of sea level height frequency are reported during one year, including maximum sea level heights of about 45 cm that may exceed the tide amplitudes. These can be related to storm surge events that occur only a few hours in a year, when water is pushed from the sea onto the land due to a temporary decrease in atmospheric pressure and wind.

We preferred to adopt this local vertical datum instead of the value of the geoid elevation for the Italian region provided by the International Service for the Geoid-ISG-ITG2009 (http://www. isgeoid.polimi.it/Geoid/Europe/Italy/ITG2009\_g.html) [31] since it is an independent and more accurate elevation datum. The International Service for the Geoid (ISG) estimates a geoid height in the Italian Geoid (ITG) 2009 for the coastline of Pyrgi at 48.319 m. This elevation corresponds to the contour line equal to zero in the Italian height reference frame. The DSM height reference frame is 0.42 m, as estimated by GPS/RTK data at a GCP located at the sea level along the coastline. Its elevation was corrected for the tidal range at the time of the surveys. The reference mean sea level estimated by the tidal analysis provided a local mean sea level with an uncertainty of ±11.8 cm.

#### **6. Vertical Land Motion (VLM) at Pyrgi**

The current rate of VLM at Pyrgi was estimated by the analysis of the available GPS data collected at the nearest GNSS stations of TOLF, belonging the INGV Rete Integrata Nazionale GPS (RING) network (DOI:10.13127/RING) and MAR8, belonging to the Topcon-NetGeo network (http://www.netgeo.it). These stations, which are located at about 6.5 km and 0.3 km of distance from the study site, respectively (Figure 1), have a robust time series that span for 2004–2019 for TOLF (15.23 years) and 2012–2019 for MAR8 (7.38 years) (Figure 6).

**Figure 6.** The vertical components (UP) of the time series of the GNSS station for (**a**) TOLF (time span 2004–2019, about 15 years) and (**b**) MAR8 (time span 2012–2019), both located near Pyrgi (see Figure 1 for location).

GPS data analysis has been carried out following the procedures already described in Reference [32] and updated in Reference [33] by adopting a three-step procedure using the GAMIT/GLOBK V10.7 [34] and QOCA software. This is part of a continental-scale GPS solution including >3000 stations [34]. The daily positions of TOLF and MAR8 have been estimated in the GPS realization of the ITRF2008 frame [35], i.e., the IGb08 reference. The position time series have been analyzed in order to estimate and correct offsets due to station equipment changes while, simultaneously, estimating annual and semi-annual periodic signals and a linear velocity term, whereas velocity uncertainties have been estimated adopting a power law + white noise stochastic model, as in Reference [36]. The results show that both sites are relatively stable in the IGb08 reference frame with a vertical velocity of −0.061 ± 0.135 mm/year for TOLF and −0.456 ± 0.344 mm/year for MAR8. We remark that uncertainties associated on the vertical velocities are about ±0.5 mm/year and are barely significant in view of unresolved questions about the GPS reference frame stability and additional factors [30].

In addition to GPS data, the tectonic stability of this region is also inferred from the low level of seismicity deduced from historical data [37] and instrumental recordings of earthquakes (www.ingv.it), which do not report the occurrence of significant events for the last 3000 years BP.

Lastly, assuming that the area will continue in the near future to have the same tectonic trend shown in the past, it is reasonable to neglect the contribution of VLM in the sea level rise projections and flooding scenarios for 2100 A.D.

#### **7. Relative Sea-Level Rise Projections and Flooding Scenarios for 2050 and 2100 A.D.**

To estimate the sea-level rise for 2050 and 2100 A.D. at Pyrgi, we referred to the regional IPCC AR5 sea-level projections discussed in the Fifth Assessment Report of the IPCC-AR5 [3], www.ipcc.ch (data available from the Integrated Climate data Center-ICDC of the University of Hamburg, http://icdc.cen.uni hamburg.de/1/daten/ocean/ar5-slr.html). These data consist of the sea-level ensemble mean values and upper/lower 90% confidence bounds of the sea level on a global grid (spatial resolution 1◦ × 1◦), obtained by adding the contributions of several geophysical sources driving long-term sea-level changes: (1) the thermosteric/dynamic contribution (from 21 CMIP5 coupled atmosphere-ocean general circulation models AOGCMs), (2) the surface mass balance and dynamic ice sheet contributions from Greenland and Antarctica, (3) the glacier and land water storage contributions, (4) the Glacial Isostatic Adjustment (GIA), and (5) the inverse barometer effect [1]. Projections, which are based on two different Representative Concentration Pathways RCP 2.6 and RCP 8.5 while providing the least and most amounts of future sea level rise, respectively, were used. The IPCC regional sea-level rate at the grid point closest to the location of the tide gauge station (Civitavecchia) was considered. By accounting for VLM from GPS data, very high-resolution DSM and regional IPCC sea level projections at the grid point corresponding to the investigated area, the first marine flooding scenarios for Pyrgi for 2050 and 2100 A.D. have been realized. To include the VLM effect in sea-level projections, we substituted the modelled GIA contribution to the IPCC rates with the GPS vertical velocities, which includes both GIA and tectonic components. Uncertainties for the sea-level estimations were calculated by combining lower and upper sea level bounds from IPCC projection and errors from GPS measurements. In any case, given the tectonic stability of the area, the VLM have a null contribution in the analysis. The relative sea-level rise in RCP2.6 and RCP8.5 scenarios at 2050 and 2100 A.D. with respect to the chosen reference epoch 2017 are shown in Figure 7, and numerical values are reported in Table 2.

**Figure 7.** Relative sea level with respect to the 2017 level as obtained from the regional IPCC sea-level projections, AR5 RCP2.6 (blue line), and RCP8.5 (red line) for a null VLM. Color bands represent the 90% confidence interval. The small-scale variations observed in the data are related to the ocean component contribution accounting for the effects of dynamic Sea Surface Height (SSH), the global thermosteric SSH anomaly, and inverse barometer effects (Church et al., 2013a, b, http://icdc.cen.unihamburg.de/). Given the vertical tectonic stability of the area, the VLM have a null contribution in the projections.


**Table 2.** Relative Sea Level Rise (RSLR, cm) above the current mean sea level at Pyrgi for 2050 and 2100 for the IPCC 2.6 and 8.5 climatic scenarios, in the mean and maximum high tide conditions. Given the vertical tectonic stability of the area, the VLM were not considered in the analysis.

The projected coastline positions for 2100 A.D., corresponding to the RCP2.6 and RCP8.5 scenarios, are obtained from the DSM for the sea levels listed in Table 2 and shown in Figure 8. The computed and the represented scenarios correspond to the local mean sea level (estimated with the uncertainty of ±11.8 cm), and are obtained by neglecting the periodical contribution due to diurnal and semidiurnal tides. To account for it, in order to estimate the maximum inundation scenarios, the time series of daily hydrometric sea levels at Civitavecchia tidal station from January 2010 to December 2018 were included in the analysis. The average half amplitude of the daily tides was estimated as high as about 30 cm (Figure 5b). Consequently, in the RCP 8.5 scenario for 2100, if we take into account both the sea-level rise (56 cm) and the mean daily tides (30 cm), we can infer a maximum water level of 86 cm. Since the highest sea levels may reach 45 cm for a few hours a year during temporary meteorological conditions when a decrease in atmospheric pressure and wind push water from the sea onto the land (Figure 5b), the maximum expected water level may reach up to 101 cm above the current level.

**Figure 8.** Projected coastline positions for 2100 A.D. for the 2.6 and 8.5 IPCC scenarios (see legend for details). Coastline positions include the uncertainty of ±11.8 cm in the reference mean sea level position as estimated from the analysis of tide gauge data. The cross sections along the shore are also shown in the figure.

For the RCP 2.6 scenario for 2100, considering a sea-level rise of 30 cm and the same mean daily tides (30 cm), a maximum sea-level rise of 60 cm can be inferred. During the rare maximum highest sea levels (45 cm) (Figure 5b), the maximum expected water level would reach 75 cm above the 2019 level.

The sea level rise scenarios for 2100 A.D. depicts the impact of the marine flooding for the coast surrounding the castle of Santa Severa, the temples of Pyrgi, and nearby beaches. The cross sections shown in Figure 9 traced along a direction that is normal to the coast and across the archaeological area of Pyrgi, which highlights the presence of a soft dune system running parallel to the shore. This system is well developed south of the castle for the whole bay of Pyrgi, while, in the north side, it is buried by the modern buildings of the Santa Severa village.

**Figure 9.** Return time (years) values of the significant wave heights (m).

Since the Etruscan time (2450 years BP), the coastline has retreated at about 80 m, at a mean rate of 3.3 cm/yr [38]. The temple of Pyrgi was originally placed at more than 120–125 m from the present-day shore and well protected from the sea by the dune system. Today, it is placed at about 40–45 m from the coastline and is undergoing water intrusion and frequent flooding, especially during storm surges (Figures 8 and 9). The proximity of the sea and the accelerated erosion of the dune system, which is highly exposed to waves especially during storm surges, are causing a continuous coastal erosion and beach retreat. Based on our sea level rise projections, a beach retreat of about 10 m is expected for 2100 A.D. (see cross sections B, C, and D in Figures 8 and 9).

This rapid retreat will accelerate the erosion process and the likely dismantling of the soft dune system with the consequent direct exposure of the temples to the sea. At the same time, the beach located along the north-west side of the castle, which is lacking a dune system and characterized by a gentle slope, is expected to retreat up to about 25 m (see cross section A), exposing the foot of the castle to severe erosion and flooding.

Lastly, Table 3 shows the expected flooded areas in 2050 and 2100 A.D. for the RCP 8.5 and RCP 2.6 climatic scenarios, together with the corresponding percentage of beach loss, which may reach up to 46% during the highest tides.


**Table 3.** IPCC scenarios, relative sea level rise projections for 2100, expected land loss (m2), and percentage of land loss with respect to the current surface for the given Relative Sea Level Rise (RSLR) projections. Given the vertical tectonic stability of the area, the VLM were considered in the analysis.

#### **8. Storm Surge Scenarios**

To evaluate the storm surge scenario for Pyrgi, a local climate wave assessment was performed using the Weibull distribution and the Equivalent Triangular Storm model originally proposed by Boccotti [39]. The return values of the significant wave heights and related peak wave periods (Table 4) have been evaluated on the parameters of the Weibull distribution [40] using data collected at the wave buoy of Ponza (wave data are freely available at http://dati.isprambiente.it/), which is located in the central Tyrrhenian sea in front of the coasts of Latium (at coordinates Latitude 40.867◦N and Longitude 12.950◦E).

**Table 4.** Return values RT (years) of the significant wave height HS (m) and the related peak wave period TP (s).


The climate wave data for ordinary and extreme storm wave conditions in terms of the return period (RT) for RT = 1 year and RT = 100 years in a sea level rise condition for 2100 in the RCP8.5 climatic scenario reported in Tables 2 and 3 were used to estimate wave setup and wave run-up during the event for both return times (Figure 10).

In the present analysis, wind set up was neglected and the precautionary assumption that the direction of wave travel does not produce refraction was considered. The expression derived by Weggel [41] and the relation of Komar and Gaughan estimated breaker depth index and breaker height index [42]. The breaker type was correlated with the surf similarity parameter, which allows the run-up estimation by the use of the predictive equations [43,44]. Based on our topographic surveys, a mean beach slope of 3.7% for the investigated coast, was assumed for the wave run-up assessment.

**Figure 10.** Projected flooded areas corresponding to the maximum run-up levels for storm surges with return times for: (**a**) 1 year (in blue at 3.37 m) and (**b**) 100 years (in blue at 5.76 m). The coastline position on 26 April 2019 is in black, which is estimated with an uncertainty of ±11.8 cm from the analysis of tide gauge data.

The results of the wave run-up analysis in ordinary and extreme storm conditions (Table 5 and Figure 11), were estimated at 3.37 m and 5.76 m for return time of 1 and 100 years, respectively.

**Table 5.** Values of maximum Rmax and medium Rmed wave run-up (m) estimated in the offshore wave conditions (HS and TP), beach slope (s m/m), breaking depth db (m), and self-similarity parameter (ξ0).


**Figure 11.** Cross sections of the coastline. See Figure 7 for cross section positions. A-A' (Beach North–Castle) north of the castle of Santa Severa; B-B', C-C' and D-D' (Beach South–Temples) across the beach and the three structures of the temple of Pyrgi. The expected sea levels for 2100 in the RCP8.5 scenario and in the maximum RCP8.5 high tide condition, are shown by the yellow and red lines, respectively. Light blue is the current mean sea level position. For storm surge scenarios for return time 1 and 100 years in sea level rise conditions, the area is almost entirely flooded with the possible exception across D-D', which remains protected by the high dune system. Projected coastlines include the uncertainty of ±11.8 cm in the mean sea level position as estimated from the analysis of tide gauge data.

#### **9. Discussion and Conclusions**

In this study, we assessed the marine flooding scenarios for 2050 and 2100 for the relevant heritage sites of Pyrgi on the basis of (a) a high-resolution DSM, (b) rates of VLM from GPS data, and (c) tidal data from the nearby tide gauge and the RCP 2.6 and 8.5 climatic scenarios released by the IPCC. The results highlight how the flooding might impact both the beach and archaeological sites with a potential local sea-level rise for the IPCC 8.5 scenario of about 56 cm for 2100 A.D. During the highest sea levels estimated by the tidal data at the nearby Civitavecchia station, the sea-level could potentially reach more than 86 cm and up to 101 cm. In addition, 2450 years ago, the shoreline was extended seaward of about 80 m more than today while, during the ancient Roman period, the sea level was at about 1.2 m below the current position, as estimated by the fish tank data [34]. Hence, a continuous sea level rise is occurring since the last millennia with an acceleration that started about 100 ± 53 years ago, at the beginning of the Industrial Era [25]. Due to the potential significant impacts on both the

coast and the heritage sites of the Temple of Pyrgi and the Santa Severa Castle with a beach retreat up to about 25 m, the expected scenario reported in this study can support adaptation plans at different time scales, which are in agreement with the Protocol on Integrated Coastal Zone Management (ICZM, https://eur-lex.europa.eu/eli/prot/2009/89/oj) in the Mediterranean.

Our results detail previous studies for the Italian [7–9,12,45–47] and the Mediterranean [10,48] regions and can contribute to raise awareness of policymakers and heritage managers toward the coastal hazard by highlighting the need to adapt actions to protect Pyrgi from marine flooding and erosion under the current conditions due to the expected sea level rise and storm surge scenarios.

We remark that the unceasing coastal erosion since the Etruscan time [38] and the retreat of the soft cliffs characterizing the Pyrgi coastline mostly occurs during high energy marine events. During storm surges, the waves approaching the coast from the Northwest, West, and Southwest marine sectors are particularly dangerous for Pyrgi. The waves approaching the coast with fetches up to about 400 km are weakly slowed down by the seafloor morphology, which results in an increase of wave energy in sea level rise conditions [49] and leads to a maximum wave run-up of 5.76 m for a return time of 100 years. The impact of extreme events may be significantly amplified by the effects of climate change. In this scenario, the intensities of the strongest future storms will exceed the strength of any in the past in the Mediterranean Sea and the oceans [50–53]. The effects of the sea-level rise in the projected scenarios include severe coastal erosion and potential enhanced damages to the heritage site of Pyrgi, which, as many other coastal areas of the Mediterranean, deserves rapid actions for their future preservation.

**Author Contributions:** Conceptualization, M.A. and A.V. Methodology, M.A. Software; validation, M.A., A.V., and M.G. Formal analysis, A.V., F.D., E.S., M.G., and G.M. Investigation, L.P., R.C., M.A., and F.D. Resources, M.A. Data curation, M.A., A.V., F.D., and M.G. Writing—original draft preparation, M.A., A.V., and F.E. Writing—review and editing, M.A. and M.G. Visualization, M.A. and L.P. Supervision, M.A. Project administration, M.A. Funding acquisition, M.A. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Acknowledgments:** The SAVEMEDCOASTS Project, funded by the EU (Agreement Number: ECHO/SUB/2016/742473/PREV16, coordinator Marco Anzidei) and Istituto Nazionale di Geofisica e Vulcanologia partially supported this study. We are thankful to Vincenzo Sepe of the INGV, for his contribution during preliminary GPS/RTK surveys.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


#### *J. Mar. Sci. Eng.* **2020**, *8*, 64


© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

## *Article* **UAV Photogrammetry and Ground Surveys as a Mapping Tool for Quickly Monitoring Shoreline and Beach Changes**

#### **Antonio Zanutta \*, Alessandro Lambertini and Luca Vittuari**

Department of Civil, Chemical, Environmental and Materials Engineering (DICAM), University of Bologna, Viale Risorgimento 2, 40136 Bologna, Italy; alessandro.lambertini@unibo.it (A.L.); luca.vittuari@unibo.it (L.V.) **\*** Correspondence: antonio.zanutta@unibo.it; Tel.: +39-051-2093111

Received: 20 December 2019; Accepted: 14 January 2020; Published: 18 January 2020

**Abstract:** The aim of this work is to evaluate UAV photogrammetric and GNSS techniques to investigate coastal zone morphological changes due to both natural and anthropogenic factors. Monitoring morphological beach change and coastline evolution trends is necessary to plan efficient maintenance work, sand refill and engineering structures to avoid coastal drift. The test area is located on the Northern Adriatic coast, a few kilometres from Ravenna (Italy). Three multi-temporal UAV surveys were performed using UAVs supported by GCPs, and Post Processed Kinematic (PPK) surveys were carried out to produce three-dimensional models to be used for comparison and validation. The statistical method based on Crossover Error Analysis was used to assess the empirical accuracy of the PPK surveys. GNSS surveys were then adopted to evaluate the accuracy of the 2019 photogrammetric DTMs. A multi-temporal analysis was carried out by gathering LiDAR dataset (2013) provided by the "Ministero dell'Ambiente e della Tutela del Territorio e del Mare" (MATTM), 1:5000 Regional Technical Cartography (CTR, 1998; DBTR 2013), and 1:5000 AGEA orthophotos (2008, 2011). The digitization of shoreline position on multi-temporal orthophotos and maps, together with DTM comparison, permitted historical coastal changes to be highlighted.

**Keywords:** coastline mapping; geomatics; SfM photogrammetry; network RTK

#### **1. Introduction**

The morphology of coastal areas is influenced by both natural phenomena and anthropic activities. The most characteristic feature of a coastal area is its shoreline, the boundary between land and a water surface.

Its location varies continuously over time due to different factors which modify the shape of the beach [1], which can create huge economic problems for the tourist industry. Coastal erosion has always existed and, over the centuries, has contributed to shaping coastal areas, through the delicate environmental balance between fluvial regulation and rising sea levels. The increase of anthropic activity along the coasts has damaged the equilibrium, accentuating the dynamics of erosive—progradation.

The coastal erosion affecting the entire North Adriatic coast puts at risks local infrastructure, environment and the tourist economy. Italian Regional Authorities, have adopted an "Integrated Coastal Zone Management Plan" involving strategies based on coastal protection structures, submerged breakwaters and, as preferred strategy for coastal protection, beach replenishment [2–4].

Over the years, research on strategies for coastal preservation has been carried out. In this context the role of geomatic techniques is fundamental for mapping and monitoring [5].

Existing shoreline mapping techniques are characterized by different accuracy, expense and training time requirements.

Topographic observations and remote sensing methods are the most commonly used techniques to detect the position of the shoreline, and monitor the shape and extension of coastal areas [6]. Point-based measurements are carried out by means of different measuring procedures with Levels, Total Stations (TS) and Global Navigation Satellite System (GNSS) instruments.

GNSS surveys (by static, fast static, kinematic, real time kinematic methods) and those carried out by means of TS (triangulation, trilateration, intersection—resection, traverses, radial sides shots methods and more generally 3D or 2D networks) together with geometric or trigonometric levelling require physical contact with the portion of land to be measured. The main limitation of these methods is the huge amount of time required and the technical difficulties involved in surveying.

The geomatic surveying methods, such as aero photogrammetric surveys [7,8], unmanned aerial vehicle (UAV) [9–15], remote sensing with satellite images [12,16], airborne light detection and ranging (LiDAR) [17–22], synthetic aperture radar (SAR) [23,24], video systems [25–29] are widespread techniques to measure wide areas, generate digital terrain models (DTMs) and maps, and to quantify coastal changes by multi-temporal comparisons.

The shoreline is constantly changing therefore all the abovementioned remote sensing techniques are fundamental because allow to acquire synchronous information that can be analysed in laboratory at a later stage. To verify the efficacy of the beach nourishment works, monitoring topo-bathimetric surveys are necessary.

The nourishment methods based on the artificial transport of sand in the areas subject to erosion, can be implemented by land, transporting the sand by truck or taking it from submarine borrowing areas located offshore and transporting the sand in situ by pipes. Given the high cost of handling the sand transport, the monitoring interventions must allow the best resolution, with the ease of use method and an excellent quality/cost ratio.

Due to the high mobility of the sand surface, for natural causes or following human intervention (e.g., for the protection of the beach touristic infrastructures during the winter), the level of significance of the comparison between two surveys of the beach surface can be quantified in ±10 cm vertically. This precision threshold can be achieved, quickly and economically by using UAVs combined with image matching photogrammetric procedures [30–33] over small-medium sized areas, making this combination an efficient tool for high resolution investigation for coastal management.

GNSS provides high accuracy coordinates of points and can be used to describe surface and volumetric changes. GNSS and photogrammetry are complementary because GNSS is fundamental for solving image orientation procedures and assessing accuracy, while photogrammetry can produce dense DSMs with radiometric content.

The aim of this work is to evaluate photogrammetric and GNSS techniques for performing 3D surveys of coastal environments in the context of coastline change studies. The techniques used have the advantage of carrying out an almost synchronous survey on medium-sized areas with the required resolution. Three multi-temporal UAV surveys were performed using UAVs. In order to allow a robust possibility of a quantitative comparison between the surfaces surveyed in subsequent epochs, the positions of the ground control points (GCPs), were surveyed, by means of a network real time kinematic (NRTK) service. In this way, in fact, the connection to the international geodetic reference system is guaranteed by a network of permanent GNSS stations. At the same time, postprocessed kinematic (PPK) surveys were carried out to produce three-dimensional models by means of interpolations of the tracks used for comparison and validation. DTMs, sections and orthophotos were obtained through photogrammetry processing with similar geometric resolution. Vectorial elaborates were numerically compared to understand and quantify the geomorphological changes while orthophotos were used to highlight variations in the coastal shoreline.

A test was carried out to evaluate the effects of the number and distribution of GCPs used to orient the blocks of images acquired from the UAV. A subset of GCPs arranged along the same route was adopted and the differences between the DTMs obtained were analysed. Furthermore, multi-temporal analysis was carried out by collecting a LiDAR dataset (CC BY-SA 3.0 IT) surveyed in

2013 provided by "Ministero dell'Ambiente e della Tutela del Territorio e del Mare" (MATTM) within the context of "Piano Straordinario di Telerilevamento Ambientale" (PST-A), 1:5000 Regional Technical Cartography (CTR, 1998) Regional Topographic Database (DBTR 2013), and 1:5000 AGEA orthophotos (2008, 2011) from the GeoER, GIS Office of Emilia Romagna Region, Italy ("Archivio Cartografico Regione Emilia—Romagna", https://geoportale.regione.emilia-romagna.it/it).

#### **2. Study Area**

The area surveyed is located between the villages of Punta Marina and Lido Adriano, on the Northern Adriatic coast 7 km away from the city of Ravenna, Italy (Figure 1). The southern jetty of Ravenna is 7 km to the North. The nearest river is Fiumi Uniti, which flows 3 km south of the area.

**Figure 1.** The survey area (highlighted in the red box), is located on the northern Adriatic coast, close to Lido Adriano, Emilia Romagna Region (Italy), 7 km from the southern jetty of Ravenna (**A**). The four cross-sections analysed (S1–S4) are identified by dashed red lines. The Regional Technical Cartography (2013) shown in black is superimposed on the AGEA Orthophoto (2011; EPSG:25833). The old shoreline (CTR 1998) is represented with black crosses (**B**).

The area studied is approximately 450 m long in the NW-SE direction and 110 m wide. It is a portion of beach especially chosen to test several geomatic techniques of surveying. Four cross-sections labelled S1–S4 from North to South (Figure 1) are defined in the study area for the analysis carried out in this work. S1–S4 are perpendicular to the average shoreline and equally spaced with a distance of 100 m between them.

The area is of considerable touristic interest, as the beach is an important resource for recreation, and also for research into the engineering of coastal zone management. From a geological point of view, the beach is characterised by fine sands covering muddy-clayey materials from alluvial deposits and ancient swamps.

The extension of the beach is a function of river and marine coastal processes influenced by anthropic activities and natural phenomena. The predominantly erosive tendency of the area is due to the limited supply of river sediments and is amplified by the rate of subsidence [26,34]. The position of the coastline is also influenced by tides and winds. As in all of the northern Adriatic, the tides are strongly asymmetrical, with a diurnal and semi-diurnal component. The characteristic tidal sea-level variations at Lido Adriano beach are 30–40 cm at neap tide and 80–90 cm at spring tide. The tidal effect combined with storm surges, and the action of winds can produce significant morphological changes to the beach [2].

The main winds are Bora and Scirocco [35]. The Bora is a cold, intense wind that comes from the NE and is more common in winter. The Scirocco is a warm wind that usually comes from the SE and is typical of the spring and autumn seasons and is one of the main factors responsible for events associated with high water levels. The beach is protected by offshore coastal protection structures that minimize erosive processes. To the North and South of the study area there are detached breakwaters perpendicular to the coast while offshore, semi-submerged barriers are parallel to the coast protecting the beach.

#### **3. Data Availability**

The current dataset of observations was acquired in three steps (2017, 2018, 2019) by means of DJI Matrice 600 (M600) and Spark UAVs (Figure 2, Table 1) in order to produce DTMs and orthophotos.

**Figure 2.** Panoramic image of the surveyed area and DJI M600 during the flight (**A**); UAVs used for the surveys: DJI Spark (**B**), detail of DJI M600 (**C**).


**Table 1.** Summary of the photogrammetric dataset characteristics.

\* Maximum flight height of a set of experimental flights on the same area.

These UAVs are multi-rotor systems:

DJI Matrice 600, 525 × 480 × 640 mm, max take-off weight 15.5 kg, with a max. payload of 5.5 kg and flight time of approximately 30 min, equipped with a DJI M600-X5 camera;

DJI Spark, 143 × 143 × 55 mm, lightened to a weight less than 300g and flight time of approximately 16 min, equipped with a DJI FC 1102 camera.

The on-board digital cameras provide useful imagery for aerial photogrammetric surveys and are georeferenced with an onboard GNSS receiver, which provides the cameras with positional accuracy of about 10 m. The three multi-temporal UAV flights (2017, 2018, 2019) were drawn up using an orthophoto of the area and the relative height of flight was selected to obtain a ground sampling distance (GSD) of about 2 cm (Table 1). Longitudinal and side image overlapping within photogrammetric blocks were fixed respectively at 80% and 70%. The direction of the flight strips was established within the flight planning software in order to minimize the number of images and the time of flight. The positions of GCPs were surveyed by means of NRTK technique using 50 cm diameter plywood targets, surveyed before each fly. These GCPs were visible in the images and were adopted to orient the blocks of images as described in the following paragraph. PPK surveys were carried out, for a direct survey of crossing profiles of the beach during the 2019 photogrammetric survey, for DTM accuracy evaluations (Figure 3).

**Figure 3.** A numbered GCP in the survey area, made with high-density heavy plywood (**A**); GCP survey by means of NRTK technique (**B**); PPK survey performed by mounting the surveying system in a rigid backpack (**C**).

The accuracies in positioning of RTK, NRTK and PPK methods, in our test conditions, are in the order of a couple of cm horizontally and 3–5 cm vertically. In order to carry out a multi-temporal analysis, DTM, maps and orthophotos were collected from the "Ministero dell'Ambiente e della Tutela del Territorio e del Mare" (MATTM), and from the GeoER, GIS Office of Emilia Romagna Region, (Italy; Table 2).



CTR—Regional Technical Map; DBTR—Regional Topographic Database; AGEA—"AGenzia per le Erogazioni in Agricoltura"; GSD—Ground Sampling Distance (m/pixel); GCP—Ground Control Point.

#### **4. Photogrammetric Analysis**

Dense point clouds and orthophotos were generated with a commercial software (Agisoft Metashape, Professional Edition 1.5.5, Agisoft LLC, St. Petersburg, Russia) which performs automatic tie point extraction and feature matching with bundle block adjustment [30]. The software is based on the structure from motion (SfM) algorithms [30–33]. In the first stage an algorithm is applied that detects points in the source photos which are stable for viewpoint and lighting variation. Secondly, a descriptor for each point based on its local neighbourhood is generated. These descriptors are used to detect correspondence between the photos. Intrinsic and extrinsic orientation parameters of the camera are than solved and refined by bundle-adjustment algorithm.

A self-calibration procedure is used in the image orientation process adopting a unique camera model for each project, evaluating the interior orientation parameters: Focal length, coordinates of principal point, lens distortion (K1, K2, K3, p1, p2), affinity and shear parameters. GNSS-assisted block orientation was applied using as *a priori* values the GNSS coordinates of the camera projection centre (PC) locations, stored in the ancillary information of the images contained in the Exchangeable image file format (Exif) file. A dense surface reconstruction is produced from the aligned images using processing methods based on pairwise depth map computation.

Finally, after a mesh calculation, texture mapping is applied and then source photos are projected onto the model (Supplementary Materials, Figures S1 and S2). GCPs were adopted to frame the bundle block adjustment and to register DTMs in the same datum and mapping projection (EPSG: 25833; Table 3). The orientation process was based on a set of GCPs uniformly distributed over the study area (Figure 4).


**Table 3.** Summary of the Orientation process and DTM characteristics.

<sup>1</sup> Ground Control Point; <sup>2</sup> Root Mean Square Error; <sup>3</sup> Digital Terrain Model.

The accuracy of the aerial triangulation process is summarized in Table 3. The last survey was identified as the most accurate by root mean square errors. The distribution of GCPs is decisive and at the same time critical in coastal contexts. Uncontrolled deformations in the derived models can be induced when GCPs in the study area are not distributed homogeneously, as may happen in long stretches of the coast. In our experience it was possible to evaluate this effect, by processing the dataset surveyed in 2019 a second time with only four GCPs aligned with each other. This processing introduced a rotation of the model as shown in Figure 4.

**Figure 4.** The vertical difference in 2019 dataset between the reference survey processed with all 9 GCPs and the test with only 4 aligned GCPs (ID: 11, 3, 5, 6). Axis of rotation highlighted in orange. The Regional Technical Cartography (2013) is shown in black.

#### **5. PPK and NRTK GNSS Surveys**

In order to allow georeferencing of the photogrammetric blocks within the UTM-ETRS89 grid system, nine GCPs were measured before the acquisition of the images, with a Topcon Hiper V Series GNSS receiver. NRTK mode using a correction mode based on virtual reference stations (VRS) through the service Topcon, Netgeo (http://www.netgeo.it/) was adopted.

Easily transportable GCPs were constructed from of 50 cm diameter plywood targets, printed with a black and white pattern and a code for easy detection. They were distributed uniformly over the area to be surveyed before flights.

The number of GCPs was chosen to minimize survey time and therefore costs. At the same time the total number of GCPs was significantly higher than the minimum necessary. This allowed a more rigorous control of the study area and redundancy in the photogrammetric orientation procedures using some of the surveyed coordinates as check points. The expected precision is 2–3 cm in horizontal coordinates and 5–10 cm for elevation, also considering the precision of the geoid model used to transform the ellipsoidal heights into orthometric ones (Geoid model ITALGEO 2005) [36].

The NRTK stationing over GCPs was fixed at two minutes of static acquisition for each point. PPK surveys were carried out to produce three-dimensional models of the topographical surface, by mounting the surveying system in rigid backpacks (GPS rover), enabling it to track the 3D profile followed by the operators.

A reference GNSS station was established close to the beach, to collect GNSS data at 1s sync rate, throughout each survey. For georeferencing of PPK surveys, the coordinates of the master station were estimated in the ETRF2000 system through the classical static differenced approach using the simultaneous acquisitions carried out at the closest EUREF reference station MSEL (Medicina, Bologna).

Then the processing of the kinematic data was performed with the open source software RTKLIB v 2.4.3 (www.rtklib.com). Considering the small extension of the studied area (a few hundred metres) absolute errors of a few centimetres characterised the survey. For the PPK tracks, the positions of the GNSS rover stations were evaluated in kinematic mode starting from the known master station coordinates, obtaining a single position for each acquisition time (1 s).

Latitudes and longitudes were then projected in the EPSG 25833, using ITALGEO 2005 geoid model to transform the ellipsoidal heights into orthometric ones. The sparse points were then interpolated using a linear Kriging approach.

#### **6. Empirical Accuracy Assessment of PPK and Photogrammetric Surveys**

In order to check the internal accuracy of the PPK surveys carried out during the 2019 survey, the crossover error analysis (CEA) was performed on the tracks surveyed, according to Hsu [37]. The procedure is based on height determination and comparison at the intersections of PPK track lines [38].

A routine was elaborated by means of QGIS software to obtain precise height values along the track intersections. The procedure used was articulated in several steps. Firstly, each 3D point of the track was connected by lines ("linestring" function) and transformed in a Well-Known Binary (WKB) format. Then, by means of line Intersection function, each line crossover was detected. For each planimetric intersecting PPK line, a height interpolation was performed at crossover point. Finally, height differences were computed at crossover points.

As result, the histogram of the observed residuals exhibits an almost normal distribution with a mean value close to zero and a RMSE of 0.03 m for 697 intersection points (Figure 5).

The dense static clouds of points recorded during operator stops in the tracks were eliminated manually in order to optimize comparison. Statistical results (Figure 5C) reveal that these data are consistent with the expected accuracy for kinematic surveys in walking mode with antenna mounted on backpacks and so are sufficiently accurate to validate photogrammetric DTMs.

Then the accuracy of the photogrammetric survey was evaluated firstly by comparing the GCP coordinates extracted manually from DTMs with those from NRTK surveys as reference. The statistical results are shown in Table 4.


**Table 4.** Statistics of the differences between the coordinates of the GCPs and those measured on the photogrammetric DTMs.

These values are affected both by mismatching errors and the local density of DTMs.

PPK GNSS surveys provided the second opportunity to evaluate the accuracy of the 2019 photogrammetric DTM. The assessment was carried out both by computing the height difference of each PPK point with respect to photogrammetric derived DTM, and calculating the difference between two DTMs: The first derived by interpolating PPK points and the second obtained by photogrammetric processing (Figure 6).

In particular, the first comparison was made to minimise the effects of interpolation, as points located along the tracks were directly compared with the corresponding ones located on photogrammetric DTM, while the second test was made by interpolating PPK points with a kriging algorithm, then transforming them in raster format with cell size of 10 cm. Statistical data and the range of differences shown in Figure 6 confirm good agreement between the models. In fact, both approaches provided similar results.

Statistical values of cases A and B (Figure 6) are similar because the area surveyed is predominantly flat, so as a result the spacing between the PPK tracks was sufficient to represent the shape of the beach.

Shorelines were digitized on orthophotos and technical regional maps (Emilia Romagna Region, CTR and DBTR 1:5000 scale) and analysed from a qualitative and quantitative point of view.

The photogrammetric model of 2019 was then compared with LiDAR DTM and the differences evaluated numerically.

**Figure 5.** Results of the Crossover Error Analysis of track intersections of the 2019 survey. (**A**) map of the 2019 PPK tracks showing the height residuals in different colours, M mean value, RMSE Root Mean Square Error; (**B**) detail of a portion of the surveyed profiles with differences at the crossing points in mm (the same colours of Figure 5A); (**C**) Histogram distribution of the height differences expressed in mm.

**Figure 6.** Map of the differences between the PPK and photogrammetric DTM. (**A**) height difference over PPK points; (**B**) NRTK assisted- photogrammetric DTMs comparison. The red box highlights the survey area. The GCP location is shown by the red dots.

#### **7. Results**

The shoreline of the studied area, located at the village of Lido Adriano in Ravenna province is subjected to natural phenomena and anthropic actions which have led to the construction of groins to protect the beach [39,40]. Both shoreline and sand surface are constantly changing. UAV combined to GNSS techniques may be an efficient tool for coastal management. Surveys carried out by means of these techniques are synchronous, characterised by the necessary accuracy and quickly to perform.

The digitization of shoreline position on multi-temporal orthophotos and maps, together with DTM comparison, permitted historical coastal changes to be highlighted (Figures 7 and 8). Distances between shorelines were determined along the S1–S4 transects using the most recent measurement as reference, shown in Table 5.

**Figure 7.** Multi-temporal shorelines: (**A**) Key map of the layers Orthophoto AGEA 2011 and DBTR 2013 showing the studied area; (**B**) Detail image of the UAV derived Orthophoto 2019 and the DBTR 2013.

**Figure 8.** Multi-temporal comparison and transect analysis in Lido Adriano area: (**A**) difference between 2019 UAV derived DTM and 2013 LIDAR DTM; (**B**) cross-shore profiles of the beach (S1–S4), 2019 UAV derived DTM (red) and 2013 LIDAR DTM (blue).

Maximum recession between shorelines is recorded comparing CTR to DBTR. AGEA 2008 and 2011 shorelines remain essentially stable in the central and southern parts while a small difference is detected towards the North. Accretion trends are observed matching these shorelines with the UAV 2019 survey. A small erosional effect can be observed in the studied area from 2017 to 2019 (Figure 7).

In Figure 8A the differences obtained by comparing the DTM detected with 2019 UAV with the 2013 LiDAR model are shown in different colours. The variations in the central part of the beach show the seasonally man-made actions carried out to flatten and enlarge the beach.

The positive values (green), clearly show the operations of beach nourishment. Not significant changes can be seen along the dune parallel to the shoreline (orange). Negative values on the breakwaters perpendicular to the shoreline can be observed in the period considered. In the 2019 model these blocks of rocks were lower, with different areal extension and show a rate of settling of a few cm per year. For these comparisons high vegetation was removed from the UAV derived model.

**Table 5.** Differences expressed in m between the multitemporal shorelines along the four cross-sections (S1–S4) and statistical data. The most recent shoreline detected with the 2019 UAV survey, was used as reference. Values are in metres. CTR—Regional Technical Map; DBTR—Regional Topographic Database, AGEA—"AGenzia per le Erogazioni in Agricoltura".


#### **8. Discussion and Conclusions**

Beach morphology and long-term shoreline changes derive from the sum of natural phenomena and anthropic activities. The coastal shoreline is the indicator traditionally used to define the trend of the sandy coasts, highlighting depositional or erosional phenomena. In beaches subject to anthropic activities like those of the Northern Adriatic coast, this interpretation becomes invalid, because shorelines migrate landward or seaward depending not only on changing sea-level or subsidence of coastal regions, but also on the presence of low-crested structures and nourishment carried out periodically for beach protection (Arpae Emilia-Romagna, www.arpae.it).

In order to use beaches for tourism it is fundamental to preserve their width. Therefore, monitoring morphological beach changes and coastline evolution trends is necessary to plan efficient maintenance work as well as replenishment and constructing engineering structures to avoid coastal drift.

This study highlights that changes in the shoreline derived from orthophoto and regional maps do not follow a constant positive or negative rate because they are due to human intervention designed to protect the coast. Since the 1980s a number of defence structures have been built and replenishment repeatedly carried out. The artificial supply of sand, and the periodical flattening of the beach by bulldozer have contributed to decreasing erosional rates and widened the beach [2,41,42].

The purpose of this study was to test fast and cheap geomatic techniques for coastline mapping and detecting shoreline changes. Three UAV photogrammetric surveys integrated with GNSS measurement of GCPs and PPK tracks on a beach protected by a system of breakwaters located in Lido Adriano, Italy was presented.

The images acquired by UAV were elaborated with image matching algorithms [30–33] that allowed maps to be produced automatically and DTM to be extracted with high spatial resolution and accuracy.

In addition to investigating the shape of the beach, quantitative and qualitative analyses of the differences were carried out using freely available multitemporal maps and DTMs from the "Ministero dell'Ambiente e della Tutela del Territorio e del Mare", GeoER, GIS Office of Emilia Romagna Region, and "AGenzia per le Erogazioni in Agricoltura".

The survey procedures adopted were cheaper and faster compared to the traditional ones. The kinematic positioning techniques PPK, RTK and NRTK allow, in undisturbed operating conditions, better accuracies than 5 cm in planimetry and height. Applying these surveying techniques in the estimation of the positions of the shutter centres of the cameras, together with a cross-strip flight, and just a single GCP positioned in the central area of the block, it is possible to obtain accuracies

comparable to those obtainable through the use of a homogeneous distribution of GCPs [43]. In case of significant offsets between the GNSS antenna and the position of the camera, the presence of a synchronized onboard strapdown inertial navigation system (SINS) mechanically connected to the camera, allows to correct the orientation up to 0.10◦–0.05◦, even if the camera is mounted within a gimbal.

Low-cost, professional, commercial UAVs can be used with good results to produce maps, detect topographical changes, and estimate variations in rate and volume. This represents a tool for coastal managers and the regional authorities to improve their "Integrated Coastal Zone Management Plans". In further developments, the use of NRTK GNSS receivers onboard of the UAVs and high-performance inertial platforms directly connected and synced to the UAV cameras, may reduce the number of GCPs in the survey of long shorelines. Nowadays long-distance surveys are limited by the flight autonomy, which is significantly longer for fixed wing UAVs and the evolving flight regulations.

**Supplementary Materials:** The following are available online at http://www.mdpi.com/2077-1312/8/1/52/s1, Figure S1: 2019 UAV derived Orthophoto, Figure S2: DTM, EPSG:25833, GSD:0.05m.

**Author Contributions:** Conceptualization of the study, L.V., A.Z.; methodology, A.L., L.V., A.Z.; software, A.L.; validation, L.V.; investigation, L.V., A.L., A.Z.; data curation, A.L., L.V., A.Z.; writing—original draft, A.L., L.V., A.Z.; writing—review and editing, A.L., L.V., A.Z.; supervision, L.V. All authors have read and agreed to the published version of the manuscript.

**Funding:** This study received no external funding.

**Acknowledgments:** We are thankful to Ing. Luca Poluzzi, Ing. Luca Tavasci for their valuable help during the fieldwork activities. Special thanks to the student Ing. Riccardo Collina for his help in the GNSS-assisted block orientation test to produce Figure 4, in the frame of his second level Civil Engineering thesis (University of Bologna). The authors would like to thank MATTM, Emilia Romagna Region and AGEA for providing LiDAR data, maps and multitemporal orthophotos. The authors thank the Guest editor Donatella Dominici, the assistant editor Zara Liu, and the four referees for their constructive comments and suggestions, which greatly improved the manuscript.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

## *Article* **Shoreline Extraction Based on an Active Connection Matrix (ACM) Image Enhancement Strategy**

**Sara Zollini 1,\*,†, Maria Alicandro 1,†, María Cuevas-González 2,†, Valerio Baiocchi 3,†, Donatella Dominici 1,† and Paolo Massimo Buscema 4,5,†**


Received: 19 November 2019; Accepted: 17 December 2019; Published: 23 December 2019 -

**Abstract:** Coastal environments are facing constant changes over time due to their dynamic nature and geological, geomorphological, hydrodynamic, biological, climatic and anthropogenic factors. For these reasons, the monitoring of these areas is crucial for the safeguarding of the cultural heritage and the populations living there. The focus of this paper is shoreline extraction by means of an experimental algorithm, called J-Net Dynamic (Semeion Research Center of Sciences of Communication, Rome, Italy). It was tested on two types of image: a very high resolution (VHR) multispectral image (WorldView-2) and a high resolution (HR) radar synthetic aperture radar (SAR) image (Sentinel-1). The extracted shorelines were compared with those manually digitized for both images independently. The results obtained with the J-Net Dynamic algorithm were also compared with common algorithms, widely used in the literature, including the WorldView water index and the Canny edge detector. The results show that the experimental algorithm is more effective than the others, as it improves shoreline extraction accuracy both in the optical and SAR images.

**Keywords:** remote sensing; satellite images; synthetic aperture radar (SAR); Sentinel-1; WorldView-2; shoreline extraction; coastline extraction; active connection matrix (ACM); J-Net Dynamic; edge detection; canny edge detector

#### **1. Introduction**

The coastal environment is an extraordinary natural resource, not only from the point of view of the cultural heritage but for hosting resources that can be measured in terms of economic assets. It is a dynamic environment, subject to continuous and constant transformation. The coastal area is indeed a highly dynamic system where erosion and deposition phenomena are influenced by various factors. According to [1], the factors responsible for changes in coastal areas may be grouped into geological and geomorphological, hydrodynamic, biological, climatic and anthropogenic factors. The coastal environment, where large percentages of the global population live, change rapidly due to its dynamic nature. For this reason, the availability of up-to-date information on its state is of great interest.

The theoretical definition of "shoreline" is merely the transition between the sea and the land [2,3], and in general, between the land surface and the surface of a water body, such as an ocean, sea or

lake. The coastal zone is, instead, that area of land and water that borders the shoreline and extends landward [4]. In [5–7] other definitions about coastal areas can be found. In theory, the concept is very simple and intuitive but its application is actually a complex task because of the temporal variability of the shoreline itself. The temporal variability develops on scales profoundly different from instantaneous to secular variations, and it depends on various factors, including wave motion, tides and winds, but also erosion and deposition. For this reason, the results obtained from surveys on the field or remote sensing are generally indicators of the actual shoreline.

According to [7], there are three main subdivisions of shoreline indicators:


In this work, we will experiment with several indicators of the third type by comparing them with each other and with those of the first type. It is important to underline, however, that those things which are extracted are indicators of the shoreline and not the shoreline itself.

Several survey methodologies are used for the detection of shorelines: from traditional topographical and GNSS (Global Navigation Satellite System) surveys to remote sensing techniques (aerial and unmanned aerial vehicle (UAV) photogrammetry, video systems and satellite optical and SAR images [8,9]), as widely discussed in another recent work [10]. In the last few years, the use of optical and SAR satellite images for automatic and semi-automatic shoreline extraction has complemented the traditional approach based on surveying methods and air-photo processing [11]. Several approaches have been proposed to outline the shoreline using remote sensing, such as edge detection, segmentation and classification approaches [6].

An edge in an image is a boundary where there is a change in some physical parameters and it is like a contour between two different regions. Many approaches have been used for edge extraction; for example, the Canny edge detector [12,13] or the snake model [14]. Segmentation methods can be divided into three categories: threshold-based, region-based and edge-based [15]. Concerning classification, they can be pixel-oriented, based on individual pixel classifications or object-oriented—which instead are able to group characteristics by aggregation in similar regions or polygons. Among the pixel-oriented classification, spectral indices have been successfully used to extract shorelines. In particular, water indices (WIs) and vegetation indices (VIs) have been developed, which are obtained by combining two or more spectral bands. Although VIs are mainly used for vegetation analysis, they have also been tested for shoreline mapping purposes [10,16]. WIs have been specifically developed for water body detection [17].

In this work, an experimental algorithm, called J-Net Dynamic [18,19], was tested both on a very high resolution (VHR) WorldView-2 optical image, and for the first time, on a high resolution (HR) Sentinel-1 SAR image. The J-Net Dynamic algorithm can be considered an edge operator, as it highlights the edges of the image and is part of the active connection matrix (ACM) system, which is a new unsupervised artificial adaptive system developed by Semeion Research Institute. The system can automatically extract features of the images—edges, tissue differentiation, etc., when they are activated by original non-linear equations. ACM activation reduces the image noise while maintaining the spatial resolution of high contrast structures. The J-Net Dynamic exploits the dynamic connections between a central pixel and its neighboring pixels. The relationships among these latter pixels with their neighbors is also taken into account, so the new image is created considering the totality of these contributions at every elaboration cycle. Every single neighbor pixel participates in the evolution of the ACM, which is finally used to create the pixel matrix with the new image. A detailed description of the algorithm can be found in Appendix A.

#### **2. Study Area and Data Set**

The Abruzzo region (Italy) coastline has a length of 125 km, of which 26 km is high coast and 99 km is sandy coast; the latter, therefore, accounts for 80% of the entire coastline and more than 50% is under erosion. The northern sector, between the Tronto River and Ortona, is characterized by low coasts; proceeding southwards, the coast is formed, up to Vasto, by mainly high coasts, but, at Vasto Marina and San Salvo the coasts are low and sandy. The average width of the Abruzzo beaches is often less than 100 m, and in some cases, it does not reach 50 m. The dominant winds are the Mistral and the Tramontana. Studies performed on wave measurements demonstrated that among the states of significant sea, that is wave heights greater than 0.5 m, the most frequent waves have heights greater than two meters and the most intense wave motions have heights between 3.5 and 6 m, characterized by an occurrence frequency of less than 5% [20].

The study of this paper was conducted on the sandy coast of Ortona (Abruzzo Region) in which three erosion phenomena have been reported on the coast, four active landslides that have effects on the coast and three landslide crags on the sea [21]. For this reason, according to the guidelines for the defense of the coast against erosion and the effects of climate change [22], one of the fundamental elements of a coastal information system is the extraction of the shoreline. The potential of the J-Net Dynamic algorithm for identifying the shoreline, in automatic or semiautomatic mode, from remote sensing, has been tested on two types of images: multispectral WorldView-2 and SAR Sentinel-1 satellite images. The Digital Globe WorldView-2 sensor, launched in October 2009, was the first VHR, 8-band, multispectral commercial satellite [23] with four standard (red, green, blue and near-infrared) and four new (coastal, yellow, red edge and near-infrared) spectral bands [24]. The images are provided with a panchromatic band (0.5 m) and eight multispectral bands (2 m) [25]. Moreover, depending on the level of processing, the images can be supplied as "Basic," "Standard," "Ortho Ready Standard Imagery" or "Stereo Pair" products [26]. A WorldView-2 image acquired on 29 June 2010 was analyzed in this work. This image was projected to UTM WGS84 reference system (WGS84/UTM Zone 33 N, EPSG: 32633). After pre-processing, described in Section 3.1, a pan-sharpening pre-treatment was also performed to merge the high geometric resolution of the panchromatic image with the multispectral bands. A Sentinel-1 SAR image has also been analyzed in this work. The Sentinel-1 mission is the European Radar Observatory for the Copernicus joint initiative of the European Commission (EC) and the European Space Agency (ESA) [27]. Sentinel-1 images are freely available and are also supplied with a level of processing that guarantees immediate use: level-0 provides raw images, level-1 Single Look Complex (SLC), level-1 Ground Range Detected (GRD), level-2 Ocean. Seven Sentinel-1 Interferometric Wide (IW) Level-1 (GRD) images (see Table 1) have been downloaded from Copernicus Open Access Hub website [28].


**Table 1.** List of the Sentinel-1 images used in this analysis.

They consist of focused SAR data that have been detected, multi-looked and projected to ground range using the Earth ellipsoid model WGS84. The ellipsoid projection of the GRD products is corrected using the terrain height specified in the product general annotation. The terrain height varies in azimuth but is constant in range (but can be different for each IW sub-swath) [29]. The IW swath acquires data with 250 km swath at 5 m by 20 m spatial resolution. IW captures three sub-swaths using Terrain Observation with Progressive Scans SAR (TOPSAR) [30]. Ground range coordinates are the slant range coordinates projected onto the ellipsoid of the Earth. Pixel values represent detected amplitude. There is no information about the phase. The resulting product has approximately square resolution pixels and square pixel spacing with reduced speckle but reduced spatial resolution. For the IW GRD products, multi-looking is performed on each burst individually. All bursts in all sub-swaths are then seamlessly merged to form a single, contiguous, ground range detected image per polarization. Seven Sentinel-1 images are processed in order to better remove the speckle, although just the one in the middle was processed for shoreline extraction.

Information regarding tides was also taken into account to compare the instantaneous shoreline [31]. In Ortona area, the maximum high tide recorded in the tide tables is 0.6 m and the minimum height is –0.2 m, referenced to mean lower low water (MLLW). As the wave action on the Mediterranean coast is limited [32] the images can be considered at the same tidal moment.

#### **3. Methodology**

#### *3.1. Shoreline Extraction with WorldView-2 Image*

The WorldView-2 images are provided as Ortho Ready Standard products, so a pre-processing is required to correct the geometric distortions [16]. An orthorectification was performed using ERDAS Imagine 2015 [33], and rational polynomial coefficients (RPCs) were used to correct the geometric distortions. In addition, 17 ground control points (GCPs) and 8 check points (CPs) were used to improve the process [10]. Finally, the original spatial resolution of the multi-spectral bands was improved with the panchromatic one by applying the HCS (hyperspherical color dpace) resolution merge algorithm, designed for WorldView-2 [34] and implemented in Erdas Imagine. After the pre-processing, two algorithms were applied to extract the shoreline edges: the WVWI (WorldView Water Index) and the J-Net Dynamic. The WVWI (1) is a WI adjusted for the WorldView-2 images [35] and defined as:

$$WWWI = \frac{CB - NIR2}{CB + NIR2} \,\prime \tag{1}$$

where *CB* and *NIR*2 are referred to the Coastal band (400–450 nm) and Near Infrared 2 band (860–900 nm), respectively. WVWI combines the Coastal and the Near-Infrared 2 bands pixel by pixel and it ranges mathematically between −1 and 1. WVWI tends to 1 in the presence of water and to 0 where sand is present. To vectorize the shoreline extracted from WVWI, a threshold was set in correspondence of the inflection point evaluated on some WVWI profiles, as it was determined in [36] and in a previous work [10] (Figure 1). The peaks, showed around 150 m on the profile, correspond to the rocks but, for the purpose of this paper, they were not taken into account.

A reclassification based on the threshold was performed and finally the shoreline was digitized with the Arcscan tools in ArcMap (ESRI, Redlands, California). Finally, a smoothing algorithm was used to improve the quality of vector line [37].

The J-Net Dynamic algorithm has been tested on a single band or band combinations. For the sake of brevity, in this work only the analysis of the NIR2 band, which gave the best results, is reported. The NIR2 band is indeed suitable for discriminating water from wet and dry sand, as it is almost completely absorbed by the water and possible perturbations are minimized by the shallow bottoms [38,39].

**Figure 1.** Spatial profiles of WVWI (WorldView Water Index): along the x-axis the distances (m) and along the y-axis the WVWI pixel values; the vertical line represents the inflection point in which the threshold between water and sand has been set.

The J-Net Dynamic algorithm, thoroughly described in the Appendix A, permits one to obtain a clear separation threshold between water and land (Figure 2).

**Figure 2.** J-Net algorithm tested on NIR2 band of the WorldView-2 image.

The vectorization of the shoreline extracted by J-Net is performed in the same way of the WVWI.

#### *3.2. Shoreline Extraction with Sentinel-1 Image*

The pre-processing of Sentinel-1 images is a key issue of SAR elaboration because, for a better identification of the shoreline, it is essential to improve the image quality, and in this case, to reduce noise while preserving the edges. Pre-processing methods can be grouped into noise reduction and image correction [6]. SAR images are affected by speckle, an effect caused by the coherent radiation used by radar systems. It is like a salt and pepper effect that happens because each resolution cell associated with an extended target contains several scattering centers whose elementary returns, by positive or negative interference, originate light or dark image brightness [40]. There are some adaptive filters, such as the median filter, which are suitable for SAR speckle removal. The median filter uses the median values within a moving kernel in place of each pixel of the image. That way, it smooths the image without smoothing the edges [6]. The SNAP (Sentinel Application Platform) software is used to pre-process the Sentinel-1 SAR images. The pre-processing steps that were carried out are:


The pre-processing aims to enhance the image for a better interpretation. In particular, the purpose is to increase the contrast while preserving the edges, using bands, algorithms and polarizations which maximize the difference between coast and water. The following options were tested:


Two masks, one on the sea and one on the mainland, have been created and statistics have been computed using the ENVI software (Harris Geospatial Solutions, Broomfield, Colorado, United States). For each study, the digital number (DN) minimum, maximum, mean and standard deviation are calculated. The results show that the contrast between sea and land is higher using the decibel bands and VH polarization. Among the filters computed with SNAP to reduce the speckle, the one which maximizes the contrast is the IDAN, but it presents the problem of excessively smoothing the edges. For this reason, the Lee filter is used, because it is a good compromise between keeping the spatial resolution and preserving the edges [41–43].

Finally, the image selected for analysis and shoreline extraction was the one in decibels, with the VH polarization and whose speckle had been reduced using the Lee filter. Figure 3 shows the image acquired on December 31, 2016 after the pre-processing.

**Figure 3.** Pre-processed and georeferenced Sentinel-1 image. Date: December 31, 2016. Units: decibels. Polarization: VH polarization. De-speckle filter: Lee.

After image pre-processing, edge detection algorithms are used to extract the shoreline. It is also manually vectorized to create a reference. As reference shoreline, both an orthorectified Sentinel-2, acquired on 1 January 2017, and the Sentinel-1 on 31 December 2016 itself were used. Due to the inherent characteristics of the SAR images, a better reference line can derive from orthoimagery. This is the reason why the Sentinel-2 optical image was also used as reference. As it was given already orthorectified, it was pre-processed by applying the atmosphere correction using the Sen2Cor tool in SNAP, and then, it was resampled to 10 m. Finally, for the same reason previously explained in Section 3.1, the NIR band (that is band 8 in Sentinel-2) was chosen to extract the reference shoreline. Then, two edge detection algorithms were tested. The first one is the Canny edge detector, which detects a wide range of edges in raster images and produces thin edges as a raster map. It operates, at first, with a Gaussian filter (based on normal distribution) to reduce the noise, then two orthogonal gradient images are computed. Finally, only relevant or significant edges are extracted by thresholding with hysteresis [44]. The result is shown in Figure 4.

**Figure 4.** Canny algorithms applied on the pre-processed Sentinel-1 image acquired on 31 December 2019.

Then, the J-Net Dynamic algorithm was tested. The result provided by the J-Net Dynamic algorithm on the SAR image is shown in Figure 5.

#### **4. Results and Discussion**

The shorelines extracted from the WorldView-2 image are presented in Figure 6: the red line illustrates the shoreline extracted manually, while the yellow one is the one generated using the WVWI and the blue one from the J-Net Dynamic.

**Figure 6.** WorldView-2 image with three shorelines: the red line is the manually extracted shoreline; the yellow one is generated using the WVWI and the blue one from the J-Net Dynamic algorithm.

The three instantaneous shorelines extracted using the Sentinel-1 image are shown in Figure 7: the red one is the shoreline extracted manually, the yellow one is from the Canny algorithm and the blue one is from the J-Net Dynamic algorithm. The additional green line is the shoreline manually extracted from Sentinel-2, acquired on 1 January 2017; that is, one day after the Sentinel-1. It was also used as a reference for the reasons explained in Section 3.2.

**Figure 7.** SAR Sentinel-1 image on the background showing four shorelines: the red and the green lines are, respectively, the shorelines manually extracted from Sentinel-1 and Sentinel-2 images; the yellow one is the one generated using the Canny algorithm and the blue one from the J-Net Dynamic algorithm.

In order to compare the accuracy of the obtained shorelines, an index *I* (2), reported in literature [16,45], is used:

$$I = \frac{S}{L'} \tag{2}$$

where *S* is the area calculated between the shoreline extracted from the algorithm and the shoreline used as reference; and *L* is the length of the reference line. The index is expressed in meters. The values obtained are shown in Table 2.


**Table 2.** I Index values. Comparisons.

The results of the index indicate that the use of J-Net Dynamic algorithm improves the shoreline extraction both in the optical and in the SAR images compared to the common filters. For the optical WorldView-2 image, the index I is equal to 1.67 m using WVWI indexed image and 1.20 m using the J-Net Dynamic derived image. This means that, in the VHR image, the J-Net Dynamic algorithm improves the accuracy by about 28%. For the SAR Sentinel-1 image, considering as a reference the shoreline extracted from the Sentinel-1 image, the index I is equal to 12.52 m using Canny derived image and 5.86 m with the application of the J-Net Dynamic filter. In brief, in the HR SAR image, the accuracy is increased by about 58%. Considering as reference the shoreline extracted from Sentinel-2 image, the index I is 66.49 m for the Canny and 60.77 m for the J-Net Dynamic, which means an increase of almost 9%. All the results proves that the shoreline extracted from the J-Net Dynamic derived image is over-fitted to the reference compared to the ones extracted from the common filters. Two more things can be also noted: the first one is that there is a higher difference between the SAR image values and the optical image values. This happens because the SAR image has a lower resolution. The second one is that, with respect to the SAR image, there is a big difference depending on the chosen reference shoreline. The ones extracted with the algorithms are further away from the Sentinel-2 than the Sentinel-1 reference shoreline. The SAR image, due to its nature, could have little backscattering from the wet sand, which means that, in this case, it "confuses" the wet sand with water. However, the purpose of this paper is to test the J-Net Dynamic on optical and SAR images and to compare the results to the common algorithms. In all of the aforementioned cases, it proved to extract a shoreline closest to the reference one.

To better determine the accuracy of the shoreline extraction, a method already used in literature [10,46] was used. Ten transects (Figure 8), about one every 100 m and drawn from the land to the sea, were used to calculate the intersection between them and the seven shorelines extracted from WorldView-2, WVWI, J-Net Dynamic on the VHR image, Sentinel-1, Sentinel-2, Canny and J-Net Dynamic on the HR image.

**Figure 8.** Ten transects every 100 m, drawn from the land to the sea, to calculate the intersection between them and seven shorelines extracted from WorldView-2, WVWI, J-Net Dynamic on the VHR image, Sentinel-1, Sentinel-2, Canny and J-Net Dynamic on the HR image.

Then, the differences between the reference shorelines and the ones extracted from the algorithms were measured. In the end, the mean values and the standard deviations of every difference were calculated (Table 3).

**Table 3.** Statistic table. The means and the standard deviations have been calculated on the difference between the reference shorelines with the algorithms we tested. The difference has been calculated on the values obtained from the intersection between the transects and the seven shorelines.


As shown in Table 3, for the VHR image, the standard deviation is 1.13 m and the mean is 1.87 considering the common filter. On the other hand, the standard deviation of the J-Net Dynamic is 1.12 m and the mean value 1.45 m. For the HR image, considering the Sentinel-2 optical image as a reference, the results show that the mean value of the difference with the Canny derived shoreline is 70.11 m and the standard deviation is 29.20 m. Compared with the J-Net values, that are 66.09 m and 19.57 m respectively, they are higher. The same happens considering the Sentinel-1 as a reference: the mean values for Canny and J-Net are 17.21 m and 10.61 m and the standard deviations are 9.10 m and 8.21 m, respectively. It can be noticed that the mean and standard deviation are lower for the VHR optical image than the HR SAR image, as expected. Moreover, the values obtained for the SAR image using two different reference shorelines confirm the better results obtained by using the Sentinel-1 image as references, as they were found previously. Indeed, the standard deviation for the difference between Sentinel-2 and Canny shorelines is 29.20 m, compared to 9.10 m obtained by the difference between the Sentinel-1 and Canny shorelines. The same happens for the differences between the Sentinel-2 with the J-Net (19.57 m) and Sentinel-1 with the J-Net (8.21 m). However, in every case, the J-Net Dynamic gives best accuracy compared to the common filters, as already found out with the index I.

In [16], the index I was used to compare the results obtained for the shoreline extraction using the NDVI and the NDWI (normalized difference water index, which considers the same bands of WVWI and the same equation) filters on the WorldView-2 initial and pansharpened image. It confirms that NDWI application provides better results; and that the pan-sharpening permits one to enhance geometric resolution and to reduce, to less than 1 m, the mean value of shifts between the automatically extracted shoreline and the manually vectorized one. They found an index value of 1.386 m for the multispectral image and 0.955 m for the pan-sharpened one. The better accuracy, compared to the one obtained for the research of this paper, that is, 1.67 m, could be due to the fact that they identified a different threshold to detect the shoreline based on a classification method using three contrasting classes (sea, land and vegetation) or due to the fact that the range of shoreline was both sandy and rocky. The rocky parts are less sensitive to instantaneous dynamic change than the sandy ones. The index *I* was also considered to evaluate the average shift between two temporal consecutive shorelines [11]. They compared shorelines during 2005 and 2011 and they evidenced that the effects of the human intervention, like the installation of breakwaters, limit erosion phenomena near shore zone. This tendency was also remarked on for the years 2011–2012.

The index *I* has not been used yet in literature to assess the SAR images. In [47], an automated method of the shoreline position detection using Sentinel-1 SAR image was studied. The method aimed at the automatic extraction of shoreline edge from pre-processed images. The algorithm performs four steps: despeckling, binarisation, morphological operations and edge detection by means of the Canny edge detector. Once the shoreline is extracted, it is compared to the one extracted from the images acquired by video system. The VDS (video-derived shoreline) is a collection of shoreline points. So, the distance between each point from the SDS (satellite-derived shoreline) is calculated in order to

evaluate quantitatively the correctness of the SDS position. The average offset turns out to be about 10 m and the root mean square difference (RMSD) is 12.48 m.

To compare the accuracy with the case study discussed in this paper, the vertices of the reference shoreline (so the shoreline manually outlined on the pre-processed Sentinel-1 image) are extracted and the shortest distance between them and the Canny derived shoreline is calculated. The results show that the average offset is around 16 m, while the standard deviation is 11.40 m. So, a robust contrast can be performed also by applying the binarisation, as proposed by [47], because it splits exactly water and land, and as consequence, the edge detector can easily extract the contour of the image. To increase the accuracy, [48] propose an integration between RASAT pansharpened image and Sentinel-1 SAR image to improve the quality of the results. The first land/water segmentation was obtained using RASAT image by means of random forest classification method. Then, the result was used as training data set to define fuzzy parameters for shoreline extraction from Sentinel-1A SAR image. The manually digitized shoreline was used as a reference and the accuracy assessment was performed by calculating perpendicular distances between reference data and extracted shoreline by the aforementioned method. The mean distance value between the final result and the reference data was calculated to be 5.59 m, which is half pixel size of Sentinel-1A and which is almost three times better then the results obtained in this research. The higher accuracy was due not only to the integration between optical and SAR images, but also to the fact that the coast considered in [48] is mainly stony rather than sandy.

Due to the high resolution difference, the authors find it inappropriate to compare the WorldView-2 optical image and the Sentinel-1 SAR image. The purpose of the paper was only to discuss potentiality of the J-Net Dynamic algorithm both in optical and in SAR images independently. So, an additional test was performed. To compare the accuracy with the average offset obtained by means of the Canny edge detector (16 m with an RMSD equal to 11.40), the vertices of the pre-processed Sentinel-1 image shoreline were extracted and the shortest distance between them and the J-Net derived shoreline was calculated. An average offset around 7 m and a standard deviation of 7.36 m were found. The accuracy was lower than one pixel. The values confirm the trend that was delineated by means of the index I value: the algorithm J-Net Dynamic increased the accuracy of the extraction.

In future studies, other images, comparable in terms of resolution, such as Sentinel-1 and Sentinel-2 or WorldView-2 and COSMO-SkyMed images, can be integrated to improve the accuracy of the results. Moreover, to answer to the question of whether the Sentinel-1 is capable of monitoring coastal changes, another interesting future work could consider all the Sentinel-1 image data set to study the shoreline changes during a determined period of time.

#### **5. Conclusions**

Monitoring coastal areas with shoreline mapping allows one to analyze erosion and deposition phenomena in order to be able to propose corrective actions. In this work, an experimental algorithm, called J-Net Dynamic, has been presented to automatically extract the shoreline. It has been tested both in multispectral and in SAR images with different resolutions. Both images were pre-processed to improve their geometric and radiometric properties in order to facilitate their subsequent analysis. They were then processed with filters and indices commonly used for this kind of research, like the WorldView water index for the optical image and the Canny edge detector for the radar image, and the J-Net Dynamic. The J-Net Dynamic algorithm was mainly developed and previously applied in medical fields and it has been tested on a VHR optical image, and for the first time, also on an HR SAR image for environmental monitoring. To validate and verify the accuracy and robustness of this approach, a comparison with other typical algorithms was performed through an index *I*, which divides the area between the shoreline obtained by the algorithm and the reference one with the length of the latter, and through the use of ten transects along the shorelines. The results showed that J-Net Dynamic improves the detection, obtaining a shoreline closest to the reference, both for VHR and HR satellite images.

**Author Contributions:** These authors contributed equally to this work. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **Abbreviations**

The following abbreviations are used in this manuscript:


#### **Appendix A**

#### *J-Net Dynamic*

This algorithm was developed by Professor Paolo Massimo Buscema, mathematician, Director of the Semeion Research Centre of Science of Communication of Rome and Full Professor Adjoint at the University of Colorado (USA).

The patent concerns active connections matrix systems (ACM), according to which each image is considered as an active matrix (network) of connected elements (pixels) that develops over time. The main idea upon which this theory is based states that each digital image stores the maximum amount of information within the pixel values and their relationships. Furthermore, it is possible to

obtain important information by analyzing the reciprocal positions occupied by pixels. For a complete presentation of ACM algorithms, see [18,49].

Any digital image is a matrix made of as many rows as the pixels that determine the width and as many columns as the pixel number related to the height. Each pixel is identified by its coordinates *<sup>i</sup>* ∈ 1, . . . , *<sup>R</sup>*, *<sup>j</sup>* ∈ 1, . . . , *<sup>C</sup>*, and its brightness *<sup>L</sup>* ∈ <sup>2</sup>*<sup>M</sup>* (e.g., in the case of 256 shades of gray, M = 8). For each pixel *uij* a set containing all the surrounding pixels *I*(*uij*) named *neighborhood* can be defined. In the ACM systems, all the pixels *uxy* ∈ *I*(*uij*) are linked to the central pixel *uij* by means of the *w*(*ij*),(*xy*) connections. The systems are classified into three orders of complexity, according to the type of evolution over time. In the first order of complexity, the values of connections are initialized once, at the beginning, and then remain fixed while pixel values *u*[*n*] *ij* evolve over time until convergence, starting from the *u*[0] *ij* value, directly derived from the image. The situation is specular in the case of second order, where the pixels' values are fixed and equal to *u*[0] *ij* , while the connections values *<sup>w</sup>*[*n*] (*ij*),(*xy*) are updated at each iteration, after initializing them to values *w*[0] (*ij*),(*xy*) ≈ 0. Finally, the third order of complexity includes models in which both the pixels and the connections change over time.

J–Net Dynamic [19,49] is a ACM system with dynamic connections and units (third order of complexity). The main prerogative of this method is to consider in its equations not only the central pixel as such, with its relative neighborhood, but also as part of the surroundings of each of the pixels around it, when they are considered in turn as the central pixel. At the beginning of the process, the units *uij* are linearly scaled into the range [−1 + *α*, 1 + *α*], where *α* is a parameter to be set up. By varying the *α* value it is possible to study images in an iterative manner. The first part of computations involves exclusively the central pixel and its neighborhood *Iij*, as shown by Equations (A2)–(A7). J–Net follows the schema: update of weights (Equations (A1)– (A5), computation of the new pixel values based on the weights (Equation (A6)), update of units (Equations (A7)–(A14), re-update of weights and so on. In this case, the update of units also involves the neighborhoods *Ixy* such that *uxy* ∈ *Iij*.

$$S\_{ij}^{[n]} = \pi \cdot (r\_{ij}^{[n]})^2 \tag{A1}$$

$$D\_{ij}^{[n]} = \sum\_{(xy) \text{s.t.} \boldsymbol{u}\_{xy} \in I\_{ij}} (\boldsymbol{u}\_{ij}^{[n]} - \boldsymbol{w}\_{(ij),(xy)}^{[n]}) \tag{A2}$$

$$J\_{ij}^{[n]} = \frac{\mathcal{e}^{D\_{ij}^{[n]}} - \mathcal{e}^{-D\_{ij}^{[n]}}}{\mathcal{e}^{D\_{ij}^{[n]}} + \mathcal{e}^{-D\_{ij}^{[n]}}} \tag{A3}$$

$$
\Delta w\_{\langle ij \rangle, \langle xy \rangle}^{[n]} = - (u\_{ij}^{[n]} - l\_{ij}^{[n]}) \cdot (-2 \cdot l\_{ij}^{[n]}) \cdot (1 - (l\_{ij}^{[n]})^2) \cdot (u\_{xy}^{[n]} - w\_{\langle ij \rangle, \langle xy \rangle}^{[n]}) \tag{A4}
$$

$$w\_{(ij),(xy)}^{[n+1]} = w\_{(ij),(xy)}^{[n]} + \Delta w\_{(ij),(xy)}^{[n]} \tag{A5}$$

$$P\_{ij}^{[n]} = Sp \cdot \frac{1}{|I\_{ij}|} \cdot \left(\sum\_{(xy)s.t.u\_{xy} \in I\_{ij}} w\_{(ij),(xy)}^{[n]}\right) + O\_{P\_{\prime}} \tag{A6}$$

where:

$$S\_P = \frac{M\_P}{M\_{\varpi} - m\_{\varpi}}$$

$$O\_P = -\frac{m\_{\varpi} \cdot M\_P}{M\_{\varpi} - m\_{\varpi}}.$$

*MP* denotes the maximum value available for pixels, *Mw* and *mw* are, respectively, the maximum and the minimum weights *w*[*n*] (*ij*),(*xy*) ∀*i*, *x* ∈ {1, . . . , *R*} and *j*, *y* ∈ {1, . . . , *C*}.

$$Out\_{ij}^{[n]} = S\_0 \cdot \frac{1}{|I\_{ij}|} \cdot \left(\sum\_{(xy)s.t.u\_{xy} \in I\_{ij}} w\_{(ij),(xy)}^{[n]}\right) + O\_0 \tag{A7}$$

*J. Mar. Sci. Eng.* **2020**, *8*, 9

where, with the usual notation:

$$S\_0 = \frac{2}{M\_{\varpi} - m\_{\varpi}}$$

$$O\_0 = -\frac{M\_{\varpi} + m\_{\varpi}}{M\_{\varpi} - m\_{\varpi}}.$$

The internal activation state of each pixel is then defined as:

$$S\_{ij}^{[n]} = |Out\_{ij}^{[n]}|.\tag{A8}$$

Therefore, the closer the weighted average of the connections of each pixel with those of its surroundings is to a neutral value, 0 or 127 depending on the encoding, the higher the value of the internal activation state of the pixel itself. With the aim of defining an update rule for units, the quantity Δ*S*[*n*] (*ij*),(*xy*) is considered (Equation (A9). Then, according to Equations (A10)–(A11), the transition from *Iij* to *Ixy* takes place.

$$
\Delta S\_{\binom{[i]}{i}, \text{(xy)}}^{[n]} = -\tanh(S\_{ij}^{[n]} - u\_{xy}^{[n]}) \tag{A9}
$$

$$\log\_{ij}^{[n]} = L \cdot u\_{ij}^{[n]} \cdot \sum\_{(xy)s.t.u\_{xy} \in l\_{ij}} \left( 1 - \left( \Delta S\_{(xy),(ij)}^{[n]} \right)^2 \right) \tag{A10}$$

$$\Psi\_{ij}^{[n]} = \sum\_{(xy)\mathfrak{s}\mathfrak{t}.\mathfrak{t}.u\_{xy} \in I\_{ij}} \tanh(\mathfrak{q}\_{xy}^{[n]}).\tag{A11}$$

The delta quantities required for correction shall be calculated in the last step. It is possible to choose for two different update laws named union (Equation (A12)) and intersection (Equation (A13)).

$$
\delta u\_{ij}^{[n]} = \varphi\_{ij}^{[n]} + \Psi\_{ij}^{[n]} \tag{A12}
$$

$$
\delta u\_{ij}^{[n]} = \varrho\_{ij}^{[n]} \cdot \psi\_{ij}^{[n]} \tag{A13}
$$

$$
\mu\_{ij}^{[n+1]} = \mu\_{ij}^{[n]} + \delta u\_{ij}^{[n]}.\tag{A14}
$$

#### **References and Note**


c 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

## *Article* **Automatic Shoreline Detection from Eight-Band VHR Satellite Imagery**

#### **Maria Alicandro 1, Valerio Baiocchi 2, Ra**ff**aella Brigante <sup>3</sup> and Fabio Radicioni 3,\***


Received: 12 November 2019; Accepted: 11 December 2019; Published: 13 December 2019 -

**Abstract:** Coastal erosion, which is naturally present in many areas of the world, can be significantly increased by factors such as the reduced transport of sediments as a result of hydraulic works carried out to minimize flooding. Erosion has a significant impact on both marine ecosystems and human activities; for this reason, several international projects have been developed to study monitoring techniques and propose operational methodologies. The increasing number of available high-resolution satellite platforms (i.e., Copernicus Sentinel) and algorithms to treat them allows the study of original approaches for the monitoring of the land in general and for the study of the coastline in particular. The present project aims to define a methodology for identifying the instantaneous shoreline, through images acquired from the WorldView 2 satellite, on eight spectral bands, with a geometric resolution of 0.5 m for the panchromatic image and 1.8 m for the multispectral one. A pixel-based classification methodology is used to identify the various types of land cover and to make combinations between the eight available bands. The experiments were carried out on a coastal area with contrasting morphologies. The eight bands in which the images are taken produce good results both in the classification process and in the combination of the bands, through the algorithms of normalized difference vegetation index (NDVI), normalized difference water index (NDWI), spectral angle mapper (SAM), and matched filtering (MF), with regard to the identification of the various soil coverings and, in particular, the separation line between dry and wet sand. In addition, the real applicability of an algorithm that extracts bathymetry in shallow water using the "coastal blue" band was tested. These data refer to the instantaneous shoreline and could be corrected in the future with morphological and tidal data of the coastal areas under study.

**Keywords:** WorldView-2; Abruzzo; multispectral classification; shoreline; coastline

#### **1. Introduction**

The coastal environment is an extraordinary natural and economic resource that is subject to continuous transformation. The coastal area is a highly dynamic system where erosion and deposition are influenced by various factors, including meteorological, biological, geological, and anthropogenic factors. The theoretical definition of coastline is merely the transition between the sea and the land [1–3]. In theory, the concept is very simple and intuitive, but its application is actually a complex task because of the temporal variability of the coastline itself. This variability develops on scales profoundly different, from instantaneous to secular variations, and it depends on various factors including wave motion, tides, winds, erosion, and deposition. It has been noted that the most significant and potentially incorrect assumption in many shoreline investigations is that the instantaneous shoreline represents

"normal" or "average" conditions [3]. For this reason, what it is important to underline that what is obtained from survey in the field or remote sensing are generally indicators of the actual coastline [3–5]. According to some authors, there are three main subdivisions of coastline indicators: characteristics visible by an operator on an aerial or remote sensing image; intersection between a tidal datum and a digital terrain model or a coastal profile; and characteristics of multispectral images identified by automatic algorithms not necessarily visible to unaided operator [3]. In this work, we assess techniques associated with the third type by comparing them with each other and with those of the first type. The purpose of this work is to assess potentialities of automatic extraction of the coastline in a specific area particularly prone to erosion. In fact, the Abruzzo region (Italy) coastline has a length of 125 km, of which 26 km is high coasts and 75 km is sandy coasts; the latter, therefore, are 80% of the entire coastline, and more than 50% are under erosion effects. The study of this paper was conducted on the sandy coast of Ortona (Abruzzo Region), in which three erosion phenomena have been reported on the coast, four active landslides have effects on the coast, and three landslide crags occur on the sea [6]. For this reason, according to the guidelines for the defense of the coast against erosion and the effects of climate change [7], one of the fundamental elements of a coastal information system is the extraction of the shoreline. A specific assessment on these areas has not yet been carried out and its results will be very useful for the continuation of the research and for the applications that can be implemented. It is important to note that what is extracted are indicators of the coastline and not the coastline itself (Figure 1).

**Figure 1.** Indicators of the coastline.

In this work, after a brief illustration of the state of the art on the techniques of extraction of the coastline, we will introduce the materials and methods used for the experiments illustrated. Subsequently, the results obtained by them will be shown in the following specific paragraph, while the conclusions and possible future developments will be discussed together in the last paragraph of the paper.

#### **2. State of the Art of Instant Shoreline Survey**

From what has been said previously, it is easy to understand that the relief of the instantaneous shoreline is a complex feature due to the intrinsic dynamism of the shoreline itself. A wide variety of methodologies have been used [7–10] that can be summarized mainly in the following:


#### **3. Materials and Methods**

#### *3.1. The WorldView-2 Satellite*

The WorldView-2 satellite, launched in October 2009, added to the already existing constellation of commercial satellites DigitalGlobe [18], already formed at the time by the satellites WorldView-1 (launched in 2007) and QuickBird (launched in 2001). Its planned mission duration was 7.25 years, but it has just turned 10 years old. WorldView-2 operates at an altitude of 770 km with an inclination of 97.2◦ for a maximum orbital period of 100 min. It is equipped with instrumentation that allows it to collect high-resolution multispectral images and stereoscopic images.

The maximum ground resolution of WorldView-2 in the panchromatic band is 46 cm, while in the multispectral band it is 1.8 m, however the distribution and use of images with a resolution greater than 50 cm in the panchromatic band and 2 meters in the multispectral band is subject to approval by the US government. The high spatial resolution allows discrimination of details, such as vehicles, shallows, and individual trees in an orchard, while the high spectral resolution provides detailed information on different areas, such as road surface quality, sea depth, and plant health. The images are taken in eight spectral bands. In addition to the four standard bands (blue, green, red, and near infrared), WorldView-2 includes four new bands at 1.8 m resolution: coastal blue, yellow, red edge, and near infrared-2 (Figure 2).

**Figure 2.** WorldView-2 spectral bands.

#### *3.2. Test Site*

The images used in this project are multispectral WorldView-2 images representing a part of the coast of the city of Ortona (Italy) (Figure 3). For this specific project, the tests were performed on a bundle of panchromatic and multispectral images acquired simultaneously in June. To be sure that the variability of the results depended only on the algorithm used in the various tests performed and not on factors due to different inclination of the image or the sun, the tests were all repeated on the same bundle of images. The images were processed entirely on the whole scene, but the results were highlighted on a single profile that was considered significant and representative of the tests on the whole scene. Ortona is a seaside town overlooking the Mediterranean Sea, more specifically, the Adriatic Sea (Figure 4); the coast is high and rocky in the southern part while it is low and sandy in the central northern part. The inner part of the city of Ortona is partly flat and partly hilly.

**Figure 3.** The test site area on the Adriatic coast (red circle).

**Figure 4.** Test site area: the yellow box is the area shown in the following figures.

In order to evaluate the possibilities of automatic extraction of the coastline (understood as an instantaneous shoreline as visible on satellite images), different combinations of algorithms on different spectra of the available images were tested and are illustrated below. The tests were developed in the ENVI 5.5 environment [19].

#### **4. Experimentation and Results**

#### *4.1. Analysis of Individual Bands and Band Combinations*

The first step in such an analysis is almost always pansharpening [20]; this operation creates a single multilayer image combining the panchromatic image of higher spatial resolution (in this case 50 cm) with the various multispectral levels (Figure 5) that have more thematic information but lower geometric resolution (in this case 1.8 m). Several algorithms have been developed to optimize pansharpening, and the comparison between them is at the center of a scientific debate [21–23].

In this experiment, Graham Smith's algorithm [24] was used because it has demonstrated to be the most reliable and accurate in previous studies [21,23]. This algorithm first calculates the weights of the single bands with the following equations:

$$B\_{wt} = \int\_{0.4}^{0.5} OT\_B(\lambda) \ast SR\_B(\lambda) \ast SR\_P(\lambda) \tag{1}$$

$$G\_{\rm wt} = \int\_{0.5}^{0.6} OT\_G(\lambda) \ast SR\_G(\lambda) \ast SR\_P(\lambda) \tag{2}$$

$$R\_{wt} = \int\_{0.6}^{0.7} OT\_R(\lambda) \* SR\_R(\lambda) \* SR\_P(\lambda) \tag{3}$$

$$NIR\_{wt} = \int\_{0.7}^{0.9} OT\_{NIR}(\lambda) \* SR\_{NIR}(\lambda) \* SR\_P(\lambda) \tag{4}$$

*J. Mar. Sci. Eng.* **2019**, *7*, 459

where *OT* is optical transmittance and *SR* is the spectral response of the various bands with their wavelengths. Once the specific weights have been calculated, the simulated band is created using the following equation:

$$RanBand \ = (B \ast B\_{\text{wt}}) + (G \ast G\_{\text{wt}}) + (R \ast R\_{\text{wt}}) + (NIR \ast NIR\_{\text{wt}}) \tag{5}$$

**Figure 5.** Some band combinations: (**a**) NIR1-G-B; (**b**) NIR1-Rededge-R; (**c**) R-G-B; (**d**) Rededge-R-Y.

Once the pansharpening has been undertaken (Figure 6), it is possible to start to perform operations on the single bands or on combinations of them.

**Figure 6.** The pansharpening process: (**a**) MultiSpectral; (**b**) Panchromatic; (**c**) PanSharpened.

Among the most used indices in the combination of the bands there is undoubtedly the vegetation index (NDVI) and the normalized difference water index (NDWI). Both these indices have been tested because they automatically highlight the instantaneous shoreline, showing a remarkable discontinuity approaching it. The NDVI is a well-known index used to study vegetation, exploiting the different response between the near-infrared band and the red band:

$$NDVI = \frac{NIR - R}{NIR + R} \tag{6}$$

In vegetation studies, it is widely used because it allows to distinguish the plants that correctly perform the chlorophyll synthesis from those that do not. In vegetation applications it provides a dimensionless value between −1 and 1 but can also be usefully used for the coastline studies because the proximity of the shoreline assumes values between –0.5 and −0.1.

The second index, the NDWI, has been studied just to distinguish the areas covered by water from those that are not; it is actually very similar to the NDVI but uses the green band instead of the red one, and its formulation is in fact the following:

$$NDWI = \frac{NIR - G}{NIR + G} \tag{7}$$

Using the NDWI, the areas covered by water are characterized by positive values of the index while vegetation and bare soil usually show negative values. Dry sand, because of its high reflectance in the green band, usually shows positive values that are close to zero (Figure 7).

ǻ**a**Ǽ ǻ**b**Ǽ

**Figure 7.** (**a**) Normalized difference vegetation index (NDVI) and (**b**) normalized difference water index (NDWI) processed images in gray scale.

#### *4.2. Automatization of Water and Vegetation Detection*

Below, we will analyze the results that can be obtained with the various algorithms available by comparing them on a specific profile that showed the characteristic aspects of a high coast and a sandy coast. On the panchromatic image, the sea/beach limits have been vectorized on the screen by an expert photo-interpreter. The purpose of this research is in fact to verify how efficient and accurate the automatic algorithms can be compared to the manual vectorization that is currently used to extract the instantaneous shoreline from satellite images. The ground surveys, already tested in previous researches [14], are certainly more accurate than those obtained from satellite images, but on sandy shores they cannot acquire the entire line of the instantaneous shoreline at the same time as it is visible on a satellite image and, therefore, are not directly comparable with them. The beach/sea and sea/detached-breakwater limits vectorized by the operator on the profile have been considered as reference for all the algorithms tested and are highlighted by the dashed lines in final comparison of the various algorithms tested. Considering that the pixel size in panchromatic is about 50 cm, the differences between the beach/sea limits vectorized by the operator with those extracted automatically equal to or less than 50 cm were not considered significant.

In this way, it will be possible to estimate the accuracy of the various algorithms, both in absolute and relative terms. What was important to observe was the possible discontinuity in the digital numbers of the images after the processing with the various algorithms; in fact, the edge-detector algorithms can easily trace border vector lines in the presence of sharp discontinuities.

The profile executed on the image processed with the NDVI algorithm shows a strong variation in slope at the boundaries of sea/beach and sea/detached-breakwater (Figure 8).

**Figure 8.** NE to SW NDVI profile path on (**a**) RGB and (**b**) NDVI processed images; (**c**) NE to SW NDVI profile, distances on the x-axis are expressed in metres. The profile trend shows sharp discontinuities at sea/land limits.

The spectral angle mapper (SAM) and matched filtering (MF) algorithms are very similar to NDVI but use all four bands (R, G, B, IR). SAM determines the spectral similarity between two spectra by calculating the angle between the spectra and treating them as vectors in a space with dimensionality equal to the number of bands. This technique, when used on calibrated reflectance data, is relatively insensitive to illumination and albedo effects [25]. MF, instead, maximizes the response of the known endmember and suppresses the response of the composite unknown background, thus matching the known signature [26] (Figure 9).

**Figure 9.** NE to SW matched filtering (MF) muddy profile path on (**a**) RGB and (**b**) MF muddy processed images; (**c**) NE to SW MF muddy profile, distances on the x-axis are expressed in metres. The profile trend does not show sharp discontinuities at sea/land limits.

The MF/SAM ratio has the characteristic of highlighting the transition from wet to dry sand that is revealed by a rapid change in the slope of the curve that represents it (Figure 10).

**Figure 10.** NE to SW MF/spectral angle mapper (SAM) profile path on (**a**) RGB and (**b**) MF/SAM processed images; (**c**) NE to SW MF/SAM profile, distances on the x-axis are expressed in metres. The profile trend shows visible discontinuities at sea/land limits.

The availability of the coastal blue band allowed us to exploit the algorithm called relative depth, developed by Stumpf and Holderied [27]. Using combinations of all bands provided a value of water depth up to about 14 m deep.

As the name of the algorithm says, the depth is relative, and in fact the values vary between 0 and 1; for bathymetric studies, therefore, control points measured to calibrate the model are required (Figures 11 and 12).

**Figure 11.** Relative depth algorithm results represented in (**a**) gray and (**b**) pseudocolor scales.

**Figure 12.** NE to SW relative depth profile path on (**a**) RGB and (**b**) relative depth processed images; (**c**) NE to SW relative depth profile, distances on the x-axis are expressed in metres. The profile trend shows almost vertical discontinuities at sea/land limits.

Comparing the various algorithms (Figure 13), it can be observed that they all show a marked discontinuity near the actual line of contact between sea and land; it is important to note that values in the *y* direction are reported for every specific algorithm so are not directly comparable, but what is important to observe is the discontinuity in the different trends. Among them, the most efficient and accurate algorithm seems to be the relative depth, because the discontinuity shown at the passage between water and land is the deepest, showing an almost vertical trend; on the other hand, it can be noted that, at the pier, the deepness algorithm shows more than one discontinuity, making it more difficult to identify the two actual sea/detached-breakwater crossings. For this reason it was decided to deepen the study of this algorithm using the coastal blue band.

**Figure 13.** Comparison of the different algorithms on the same profile, black dashed lines are the actual water/ground limits (distances in the *x* direction are in meters along the profile, values in the *y* direction are reported for every specific algorithm so are not directly comparable).

#### *4.3. The Coastal Blue Band and the Relative Depth Algorithm*

The results obtained using the relative depth algorithm have suggested that it could be further enhanced to deepen the study by using the coastal blue band, although this is not present on all satellite platforms. One can get useful information not only on the coast but also on the first meters of the seabed thanks to its ability to penetrate clear water. For the study of the coastline, we have observed how the algorithm assumes a value of zero in the presence of land emersion (Figure 14) while reconstructing the profile of the seabed moving to the sea.

**Figure 14.** Detail of relative depth algorithm results using coastal blue band; on (**a**) the profile trac, on (**b**) the resulting profile were open water is on the right, distances in the x direction are in meters along the profile.

However, the depth is relative, and therefore some known depths are necessary to transform the relative depths into absolute values. The only official data on the depths in this area are those released by the Military Geographic Institute (IGMI), collected with various techniques and accuracy. The correct calibration of the results of the algorithms would require more accurate measurements taken at the same time as the acquisition of the image. The data provided by the IGM only allow us to make a rough evaluation of the trends of the algorithm itself. In this way, we obtained the absolute depth profile (Figure 15).

**Figure 15.** "Absolute depth" profile: along the *y* direction is the absolute depth and along the *x* direction is the progressive distance; all distances are in meters and open water is on the right.

The general trend is very similar to the expected one (Figure 16), also highlighting the differences in altitude near the detached breakwater. Despite the use of a median filter, some circular bathymetric curves are observed that might not respond to real morphologies. This could be due to the high detail of the information or to some bias of the algorithm; to be able to evaluate with certainty the performance of the algorithm, it would also be necessary to have a survey with greater accuracy (for example, side scan sonar) performed at the same time as the acquisition of the image.

**Figure 16.** "Absolute depth" bathymetric lines; open water is on the right.

#### *4.4. Supervised Multispectral Classification*

To complete the framework of the experiments, a supervised multispectral classification was also performed using the maximum likelihood algorithm; for this purpose, several regions of interest (ROI) were defined, including four building classes and three vegetation classes (Figure 17).

**Figure 17.** Detail of a vegetated area on RGB image on (**a**) and corresponding classification classes on (**b**).

The result obtained from such a classification was a new raster file in which, in each pixel, the digital number corresponding to the class of land cover to which the pixel itself has been associated is reported (Figure 18). This classification works equally well in the division of the various types of buildings and vegetation as well as with the separation of dry sand and wet sand. We compared the instantaneous shoreline thus obtained with the only official one available, that is the one reported in the regional map scale 1:5000, aware of the shifts due to different period, time, and phase of the waves, however, finding a good agreement that can only confirm that the coastline in this area has not undergone major changes.

**Figure 18.** (**a**) Pansharpened image; (**b**) classified image; (**c**) classes legend.

#### **5. Conclusions and Further Developments**

The purpose of this project was to identify the most efficient procedure for automatically extracting the instantaneous shoreline in sandy coasts of the Abruzzo region that are affected by severe erosion problems. What was important to observe were the possible discontinuities in the digital numbers of the images after treatment with the various algorithms; in fact, once any discontinuities have been identified, the edge-detector functions are able to easily draw border vector lines on the same discontinuities. Therefore, different algorithms have been tested and their trends on the sea/beach passage have been examined, comparing them with the limits identified on the same image by an experienced photogrammetric operator.

The algorithm that showed a trend with more sharp discontinuities and coinciding with those identified by the operator was the relative depth. However, this algorithm showed more discontinuities than those observable on the image at the pier. For this reason, it was subsequently experimented with to see if it was possible to improve the results using the coastal blue band. It has been verified that the costal blue with the relative depth eliminates the problem of the multiple detections. It can be deduced that the most accurate results for this specific application are obtained with the relative depth algorithm using the coastal blue band of the WorldView 2 satellite.

As far as the use of the same algorithm for the evaluation of bathymetry is concerned, it is necessary to have precise, punctual measurements to calibrate it and measurements such as side scan sonar to verify it; both should be carried out at the same time as the acquisition of the image to be significant.

**Author Contributions:** Supervision, M.A., V.B., R.B. and F.R.; Methodology, M.A., V.B., R.B. and F.R.; Investigation, M.A., V.B., R.B. and F.R.; Writing—Review and editing, V.B., R.B. and F.R.

**Funding:** This research received no external funding.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


© 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

*Article*
