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

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.


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 erosiveprogradation.
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][3][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]. Pointbased 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 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 topobathimetric 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][31][32][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, multitemporal 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).

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. 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 crosssections 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.

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.  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). 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).

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][31][32][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). 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.

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.

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.
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.

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. Table 5. Differences expressed in m between the multitemporal shorelines along the four crosssections (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". 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).

Source
In Figure

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][31][32][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 highperformance 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 www.mdpi.com/xxx/s1, Figure S1: 2019 UAV derived Orthophoto, Figure S2  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.