Mapping inland water bathymetry with Ground Penetrating Radar (GPR) on board Unmanned Aerial Systems (UASs)

Bathymetry of inland water bodies is essential for river maintenance and flood risk management. Traditionally, in shallow water bodies, bathymetry is retrieved by operators wading through the water body with Real Time Kinematic (RTK) Global Navigation Satellite System (GNSS), whilst in deeper waters, it is retrieved with sonar instruments on manned or unmanned boats. In the past, researchers have documented the use of Ground Penetrating Radar (GPR) on boats (i.e. water-coupled GPR) for monitoring the bathymetry of frozen and non-frozen water bodies. Furthermore, GPR has been used on helicopters for monitoring ice and snow thickness. However, deployment of GPR on board Unmanned Aerial Systems (UASs) in non-frozen inland water bodies with electric conductivity higher than 100 μ S/cm (as is common in most inland waterbodies in non-polar regions) is unexplored. In this paper, we document the possibility to use drone-borne and water-coupled GPR in several cross-sections located in three different waterbodies (1 lake and 2 rivers) in Denmark. These waterbodies had different bed sediment materials and vegetation conditions, an electric conductivity varying from 200 to 340 μ S/ cm and depths up to 2.5 m. Drone-borne GPR showed accuracy similar to water-coupled GPR when compared to RTK GNSS ground-truth measurements, with a Mean Absolute Error (MAE) of approx. 8 cm. The only limitations of drone-borne GPR were i) more restrictive minimum depth requirement (typically 0.8 – 1.1 m for drone-borne GPR, while 0.3 – 0.4 m for water-coupled GPR) ii) requirement to fly the GPR antenna at altitudes of approx. 0.5 m above the water surface to avoid high spreading losses and strong surface clutter events hiding the signal. Finally, GPR measurements were benchmarked against traditional sonar measurements, showing that GPR measurements significantly outperform sonar measurements in waterbodies with medium or high density of aquatic vegetation.


Introduction
Bathymetric surveys are typically conducted to measure water depth and map the underwater features of a water body.Bathymetry of inland waterbodies plays a critical role in many hydrological and hydraulic problems and applications; indeed the hydrological regime of a lake is largely dependent on its bathymetry (Yao et al., 2018) and, similarly, the predictive skills of river hydrodynamic models depend on the accuracy of river cross-section shape datasets (Grimaldi et al., 2018).
To optimize river maintenance and improve flood forecasting, bathymetry should be retrieved at regular time intervals with high spatial resolution (Alsdorf et al., 2007;Cook and Merwade, 2009).Indeed, the temporal variability of bathymetry, which is typically caused by processes such as sedimentation, erosion and human intervention (e.g.bottom dredging), can significantly change the frequency of overbank flooding under the same peak discharge conditions (Lane et al., 2007;Stover and Montgomery, 2001).Similarly, the high spatial variability of bathymetry along the river course may require continuous or very highresolution mapping to accurately characterise the river morphology and geology (Diaconu et al., 2019;Merwade, 2009).To acquire bathymetric data, different techniques can be used, as described in the next sections.

RTK GNSS surveys
In shallow water, the conventional way of measuring river bathymetry is through cross-sectional surveys with RTK GNSS at certain locations along the river; however, these surveys typically require significant man-power and are time-consuming, thus the ratio cost/area covered is very high and investigation is typically limited to few crosssections (Bangen et al., 2014;Feurer et al., 2008).

Sonar
In navigable rivers, river bathymetry is typically measured with single beam or multi-beam sonar on-board manned or unmanned vessels (Bio et al., 2020;Halmai et al., 2020;Leyland et al., 2017;Specht et al., 2020;Stateczny et al., 2019;Young et al., 2017).However, sonar systems have limitations, especially for the most commonly used frequencies (less than 1 MHz), in measuring very shallow depths due to surface clutter and multipath effects (Albright Blomberg et al., 2013).Furthermore, the accuracy of sonar significantly degrades in vegetated rivers (Helminen et al., 2019), indeed the high level of reflection of sound waves from the vegetation can result in depth measurements within vegetation canopy (Sabol, 2002).Additionally, deployment of boats can be complicated in remote areas and is limited to navigable water or, especially for unmanned vessels, locations without dense floating aquatic vegetation.To alleviate the issues of sonar deployment in remote areas or non-navigable rivers, researchers (Alvarez et al., 2018;Bandini et al., 2018) and recently companies (e.g.UgGS-SPH Engineering (Latvia) and Thurn group (UK)) developed sonar systems tethered to an Unmanned Aerial System (UAS); however, the tethered sonar has to remain in contact with water throughout the survey, which complicates automatic pilot flights and reduces the possibility to perform beyond visual line of sight flights.

Remote sensing methods
To overcome the limitations of in-situ surveys, researchers have developed remote sensing methods to estimate water depth.These methods can be categorized into indirect methods, which determine the flow area or the mean water depth, and direct methods, which determine the full bathymetric shape.Indirect remote sensing methods are based, for instance, on estimating mean depth from satellite measurements of Water Surface Elevation (WSE) and modelled discharge (Leon et al., 2006), or estimating flow area from WSE and cross-section mean flow velocity (Moramarco et al., 2019).The focus of this paper is on direct methods, which are able to directly observe submerged topography and reconstruct the full cross-section shape.Direct methods to observe bathymetry include (i) through-water photogrammetry, ii) spectral methods, iii) bathymetric LiDAR.
Spectral methods are mostly based on the wavelength-dependent exponential attenuation of light in the water column and exploit the correlation between depth and reflectance, as according to Lyzenga (1978) or Stumpf et al. (2003).Spectral methods to estimate water depth of inland waterbodies have been applied to multispectral (or hyperspectral) and RGB images from i) satellites (Geyman and Maloof, 2019;Jagalingam et al., 2015;Lyons et al., 2011;Muzirafuti et al., 2020;Wei et al., 2020), ii) aircrafts (Carbonneau et al., 2006;Flener et al., 2012;Legleiter, 2012;Marcus et al., 2003) and iii) UASs (Flener et al., 2013;Lejot et al., 2007;Rossi et al., 2020).To combine the advantages of the spectral approach, which performs better when the sediment is comparatively homogeneous (Legleiter et al., 2009;Overstreet and Legleiter, 2017), and the photogrammetry approach, which performs better if the bottom is sufficiently textured to enable feature matching and ensures higher spatial resolution (Feurer et al., 2008), hybrid approaches have been demonstrated (Slocum et al., 2020;Starek and Giessel, 2017).However, both spectral and photogrammetry methods are limited to waterbodies with high water clarity and can deliver results only for depths lower than the Secchi depth.Thus, these passive remote sensing methods are typically applied to i) shallow coastal areas or ii) in rivers with sandy or gravel bottom and water depth less than 1-1.5 m (Brasington et al., 2003;Marcus et al., 2002;Winterbottom and Gilvear, 1997;Woodget et al., 2015).
To overcome the limitations of passive sensors, active sensors such as bathymetric LiDAR systems (typically using green wavelength) have been used to monitor bathymetry of inland waterbodies from aircrafts (Bailly et al., 2010;Charlton et al., 2003;Hilldale and Raff, 2008;Kinzel et al., 2013Kinzel et al., , 2007;;Lague and Feldmann, 2020;Pan et al., 2015) and recently also from satellite missions such as ICEsat-2 (Ma et al., 2020;Parrish et al., 2019).The cost, weight and size of bathymetric LiDAR systems challenge UAS deployment.However, in the last 5 years, green LiDAR systems, which comply with the payload weight restrictions of UAS, have entered the market, such as the i) laser profiler Riegl BDF-1 (Mandlburger et al., 2016), ii) the scanning LiDAR systems ASTRALiTe Edge Lidar (Kinzel and Legleiter, 2019), iii) the scanning LiDAR TDOT green (Amuse Oneself Inc, Japan).Furthermore, Chen et al. (2021) reported the design and performance of a photon-counting bathymetric Lidar based on a fiber-laser-pumped green laser.The main limitations of the commercial LiDAR systems are i) their cost (typically more than 150 000 euros), ii) their penetration capability only up to 1-1.5 times the Secchi depth, iii) the uncertain performance for low-reflectivity bottom types and dense aquatic vegetation.

GPR
In frozen water bodies, Ground Penetrating Radar (GPR), or more accurately Water Penetrating Radar, outperforms other bathymetric remote sensing techniques because of the high penetration of GPR waves in ice and snow due to the low electric conductivity and low dielectric permittivity of these two media.In these applications, GPR can be deployed directly on the ice surface (Lapazaran et al., 2016;Liu et al., 2014Liu et al., , 2013;;Moorman and Michel, 1997;Saintenoy et al., 2013;Shean and Marchant, 2010;Yoshikawa et al., 2006) or also from helicopters (Grab et al., 2018;Kovacs, 1979;Machguth et al., 2006;Merz et al., 2015;Rutishauser et al., 2016).However, in freshwater waterbodies, the high electric conductivity and high relative dielectric permittivity of water (≈80), make GPR-based bathymetry monitoring of freshwater significantly more complex than monitoring of ice or snow.While the use of GPR deployed water-coupled from boats has been documented to successfully monitor bathymetry in lakes (Haeni, 1996;Kidmose et al., 2013Kidmose et al., , 2011;;Lachhab et al., 2015;Sambuelli and Bava, 2012;Swain, 2018) and rivers (Sambuelli et al., 2009), air-launched GPR antennas are less commonly described in the literature for bathymetric measurements of non-frozen waterbodies.Indeed, GPR penetration in deep water is dramatically reduced when the antenna is not water-coupled because of transmission, reflection and spreading losses (Delaney et al., 1990), and because the rough water surface severely degrades the coherence of the water surface transmission (Arcone and Calkins, 1990).Some tests have been performed with GPR antennas mounted on suspended cableways (Costa et al., 2006(Costa et al., , 2000;;Mason, 2004;Nobes et al., 2018;Spicer et al., 1997).However, at present, only one reference to helicopter-borne GPR monitoring of freshwater bathymetry is found (Melcher et al., 2002), although the characteristics of the sites (such as electric conductivity) were not given.Helicopter-borne deployment significantly reduces the survey time compared to water-coupled GPR; indeed, Melcher et al. (2002) reported that the helicopter could be flown at a low altitude (few meters above the surface) covering three sites located 100 km from each other within 60 min.GPR mounted on UASs would also increase survey speed compared to water-coupled GPRs and would allow for low altitude flights, which are significantly less dangerous with UASs than with manned helicopters.
In this paper, we show water-coupled and drone-borne GPR bathymetric measurements of freshwater bodies, with electric conductivity F. Bandini et al. between 200 and 340 µS/cm.We also demonstrate the advantage of using GPR compared to sonar in highly vegetated shallow (less than 2.5 m) waterbodies.

Methodology
In this chapter, we first summarize the main equations that describe the theoretical background of GPR detection of bathymetry, then we describe the survey methodology and the study areas.

Theoretical Background -GPR detection of bathymetry
The GPR signal is significantly attenuated when travelling from air into the water body, through the water body, and finally reflected from the bottom of the water body.The total loss coefficient α total estimates the ratio between the emitted power (P emitted ) and received power (P received ), or the ratio between emitted electric field-strength (E emitted ) and received electric field-strength (E received ), as shown in Eq. (1).The total loss coefficient α total describes the losses involved in GPR signal propagation from the GPR transmitter to the receiver (through the water body) and can be calculated by adding up the individual losses, as shown in Eq. ( 2), in which we have i) the transmission loss (α trans ) through the air-water interface; ii) the reflection loss (α r ) from the bottom of the water body; iii) the propagation loss (α p ) in water; iv) the spreading loss (α spread ) in air and water.In α total , we did not include losses such as i) radar system internal losses ii) temperature related losses.α total = 10log 10 (

Transmission loss
As the GPR wave propagates into water, it suffers a transmission loss (α trans ) at the surface, which can be calculated with Eq. (3), as according to Jiang and Georgakopoulos (2011), and is equal to 8.9 dB for air-water transmission.
T = 2η water η water + η air (4) In Eq. (3), multiplication by 2 is done because the wave needs to travel through the water surface twice.In the equation, η is the intrinsic impedance (computed from the dielectric permittivity), specifically η air of air or η water of water, Re{} indicates the real part of the complex number, T is the transmission coefficient given by Eq. ( 4).In the formulation of the transmission loss coefficient, we do not consider that the transmission loss may be increased because of the water roughness induced by the downwash generated by a low-altitude flying UAS; however, in real-world scenario, this roughness would indeed increase surface cluttering (Rappaport et al., 2003), affect wave coherence (Arcone and Calkins, 1990) and increase the spreading of the GPR waves exiting water, thus decreasing the amount of radiation that the receiver acquires.

Propagation loss
The attenuation of the GPR wave in water can be described through the coefficient α atten (Np/m).Eq. ( 5), which is reported in Annan (2009), shows that the coefficient of attenuation (α atten ) is a function of i) the absolute magnetic permeability of water (µ abs,water ), ii) the real part of the absolute permittivity value of water (ε abs,water ), iii) the effective electric conductivity of water (σ e water ).The absolute permittivity of water is obtained by multiplying the real component of the relative permittivity (ε water ) by the permittivity of free space (ε 0 ), similarly µ abs,water is obtained by multiplying the relative magnetic permeability of water (µ water ≈1) by magnetic permeability of free space (µ 0 ).Effective conductivity (σ e water ) can be computed through Eq. ( 6) (Cassidy, 2009), which contains the static electric conductivity of water (σ water ), the imaginary component of the absolute permittivity of water (ε ′′ abs,water ), and the angular frequency of the GPR (ω).The imaginary component of the absolute permittivity of water is obtained by multiplying the imaginary component of the relative permittivity by the permittivity of free space (Giannoukos et al., 2017).Different approaches exist to compute the GPR frequency (and thus the angular frequency) in the medium, such as using the centre frequency of the emitted wave (Lambot et al., 2004;Schwamborn et al., 2002), the entire bandwidth (Huisman et al., 2003), the frequency of the received wave or the dominant frequency expressed as 80 % of the nominal centre frequency (Grimm et al., 2006;Pérez-Gracia et al., 2008).In this research, the last approach was adopted: the GPR frequency in water was considered equal to 80 % of the nominal centre frequency and thus the angular frequency was computed accordingly.
As according to Jiang and Georgakopoulos (2011), the propagation loss (α p ) in dB/m can be calculated with Eq. ( 7), in which d is water depth and the factor two is due to the double pass in water.

Reflection loss
Reflection loss (α r ) is calculated the same way as transmission losses, but by using the reflection coefficient (R), as shown in Eq. ( 8).The reflection loss coefficient is highly dependent on the impedance of the riverbed sediment material reflecting the wave, as shown in Eq. ( 9) (Annan, 2009;Reynolds, 1997).

Spreading loss
Spreading loss, which is related to the flight height and water depth, can be calculated with Eq. ( 10), where A footprint,directWave is the footprint of the direct wave as received by the receiver, and A footprint,riverbed is the footprint at the river bottom.Computation of the GPR footprint is complicated and has been approximated by different authors with different equations (Pérez-Gracia et al., 2008).In this case, we applied Eq. ( 11) to compute the first Fresnel zone and estimate the radius at the water surface (r footprint,watersurface ) from the wavelength in air (λ air ) and the flight height (h) above the water surface, as according to Leucci and Negri (2006) and Utsi (2004).This equation is chosen because Pérez-Gracia et al. (2008) proved with their experiments in water tanks that it can give a footprint approximation better than other common equations typically used to estimate the footprint, such as the equations found in Leckebusch and Peikert (2001) or Conyers (2004), which provide a conservative estimation of the antenna's horizontal resolution capability (Pérez-Gracia et al., 2008).Experimental results determining the footprint when the GPR wave crosses multiple media (such as air and water) were not found.However, we here assume that the same equation Eq. ( 11) can be applied and thus the radius at the bottom of the water body (r footprint,riverbed ) is computed with Eq. ( 12), which is basically a sum of F. Bandini et al. the radius at the water surface after the GPR wave travels through air and the radius at the water bottom if the GPR was water-coupled.Thus, the second part of Eq. ( 12) contains also the wavelength in water (λ water ) and water depth (d).With the assumption that the footprint is a perfect circle, the area A footprint,riverbed can be computed with Eq. ( 13).A critical factor in footprint calculation is the choice of the frequency from which to estimate the wavelength in air or water.Here, for computing λ air we use the centre frequency of the Gekko-80 (≈80 Mhz), whilst, for computing λ water , we use the dominant frequency in water, i.e. 80 % of the centre frequency.

Total loss
The total loss can be computed with Eq. ( 2).An example of computation of total loss is shown in Appendix A, with a detailed calculation reported in the Supplementary Material.The computed total loss of 41.2 dB means that, out of a transmitted wave with electric field-

Table 1
Comparison of different GPR that are compatible with drone-borne deployment and suitable for bathymetry (also in non-frozen water bodies).strength of 1, only an electric field-strength of approx.0.01 can be detected in the radargram at the riverbed (i.e.rather small signal).These values do not include i) radar internal noise ii) temperature effects iii) effects of water surface roughness.Fig. 1 shows the different loss components in GPR transmission (through the water surface), propagation (through water), spreading (through both air and water), and reflection (from the bed of the water body).
In case the antenna was water-coupled, transmission loss would be approx.0 and spreading loss would be only approx.3.87 dB, indeed, a significant amount of the spreading loss occurs in air.Thus, a watercoupled GPR antenna would exhibit a smaller total loss than an elevated antenna.Specifically, a water-coupled antenna would have a total loss of 28.85 dB instead of 41.2 dB, resulting in a significant increase (12.4 dB) in return signal when the antenna is water-coupled.This is demonstrated in panels (b) and (c) in Fig. 1, in which we show the experimental result of lifting a handheld 80 MHz GPR antenna from the water surface to 1 m altitude in Furesø lake.The radargram shows the return electric field-strength of the Gekko-80 antenna when watercoupled (around trace 20) and when elevated up to 1 m (trace 430).The electric field-strength of the return wave corresponding to the lakebed (approximately 1 m water depth) decreases from 5290 (watercoupled) to approx.1190 bits when the antenna is elevated to 1 m above the water surface: corresponding to 20⋅log 10 (5290/1190) ≈ 12.9 dB loss when the antenna is elevated, thus comparable (given the uncertainty in spreading loss) to the theoretically computed value of 12.4 dB.The total loss, consisting of transmission, spreading, propagation and reflection cannot be computed (e.g. by comparing water surface and bottom electric field-strength) from Fig. 1, because manufacturer hardware gains are applied to each trace with return times greater than approx.30 ns, to improve the GPR wave penetration through water.

GPR sensors
In recent years, GPR manufactures have developed lightweight GPR units complying with the size and weight limitations of UAS and with frequencies (40-120 MHz) suitable for bathymetry measurements of non-frozen waterbodies, as shown in Table 1.
The GPR units from Radarteam and Geoscanners use a signal sampling technique called Real-Time Sampling (RTS); RTS significantly improves signal-to-noise ratio compared to the alternative interleaved sequential sampling, whilst the Malå systems use High Dynamic Range (HDR), which is similar to Real-Time Sampling.All the GPR antennas of Table 1 allow for 32-bit resolution, which results in a theoretical dynamic range of 190 dB and an effective dynamic range (after averaging, denoise, filtering) of approx.150 dB.The frequency of the antenna is a critical choice.Indeed, a lower frequency gives lower attenuation coefficient (thus lower propagation loss in water) and lower spreading loss; however, it results in worse vertical resolution (typically approximately one-fourth of the wavelength in water (Jol, 1995), thus approx.10 cm for 80 MHz) and worse horizontal resolution.

Field instrumentation
For this specific research, the GPR antenna Gekko-80 with RTS1600 data processing unit was deployed.
As sonar, the system Deeper Smart Sonar CHIRP+ (Deeper, UAB Lithuania) was deployed.This sonar records three different frequencies (each of which correspond to a different scanning angle), specifically 100 kHz (47 • ), 240 kHz (20 • ) and 675 kHz (7 • ).The sonar has a maximum depth capability of 100 m, and a minimum depth of 0.15 m (675 kHz), 0.6 m (240 kHz) and 0.8 m (100 kHz).Measurements were acquired with all the three frequencies.Although the higher frequency ensures better vertical resolution and typically higher accuracy in shallow water, the manufacturer suggests that the 240 kHz frequency is less sensitive than the 675 kHz frequency to all external factors, such as waves, wind, sonar reeling and trolling.The 100 kHz frequency was tested because it gives better penetration, which could potentially ensure better capability of penetrating aquatic vegetation.The sonar, when connected to a smartphone through the App Fish Deeper (Deeper, UAB Lithuania), displays the full waveform data in a colormap plot.However, from the logged data, the user can extract only one numeric water depth value per scan (corresponding to the value automatically estimated by the sonar software in real-time interpretation of the full waveform plot).The sonar can also be configured to send NMEA 0183 messages, containing depth and position data from the in-built GNSS receiver, via WIFI to external devices.
In addition to GPR and sonar, the deployed system included i) a GNSS receiver ZED-F9P (U-blox, Switzerland) equipped with an Antcom (3G0XX16A4-XT-1-4-Cert) multiband antenna, ii) a single-board computer BeagleBone Black (BeagleBoard.org).The single-board computer was programmed to i) log sonar data, ii) log raw GNSS observation and navigation messages and iii) send NMEA GPGGA sentences to the Gekko-80 processing unit.
The GPR system was tested both i) in water-coupled deployment on board an inflatable boat (made of non-conductive material) and ii) in drone-borne deployment on board the hexacopter DJI matrice 600 (DJI, China).The drone was equipped with a radar altimeter from UgCS (UgCS, SPH Engineering, Latvia).This radar altimeter enables flight at constant and low altitude in automatic flight missions planned with UgCS Ground Station Software.
Fig. 2 shows the water-coupled and drone-borne deployment of the Gekko-80 system.
The inflatable boat was dragged along the cross-section from one side of the river to the other (or to a specific point in the lake case study).The UAS was flown over the same cross-section while keeping a constant low altitude with the GPR antenna approx.0.5-0.6 m above the water surface.
An RTK GNSS rover station, consisting of a Trimble RTK GNSS R8s (Trimble Inc., USA) mounted on a carbon fibre pole, was deployed to

GNSS processing
The GNSS data were processed in Post-Processed Kinematic (PPK) mode with the software RTKLIB Demo5 b34b from https://rtkexplorer. com/.This specific version of RTKLIB, which differs from the standard version available on https://www.rtklib.com/(Takasu and Yasuda, 2009), enables processing of observations retrieved on the L2C GPS signal retrieved by the u-blox ZED-F9P.
The used GNSS reference station was either from a Continuously Operated Reference Station (CORS) provided by GPSnet.dk or a locally configured station (comprised of a NovAtel Flexpack6 receiver and a NovAtel GPS-703-GGG pinwheel triple frequency).

GPR processing
The GPR data were processed in the software GPRSoft™ (Geoscanners AB, Sweden), according to the following sequence of steps: • The return from the water surface is identified to define the time-zero reference point.Detection of water surface is performed with the Fig. 3. Site overview map (coordinates in ETRS89-UTM32).The maps show the river cross-sections that were surveyed in Guden Å, Nørreå, and the profile in Furesø lake.
F. Bandini et al.
GPRSoft "automated phase detection" method in which the software tries to follow a layer throughout the radargram based on the phase of a user-defined initial point, which was selected on the water surface.When the GPR is water-coupled, the water surface return is nearly a flat line.When the GPR is drone-borne, the water surface return may wiggle, because drone altitude above the water surface can fluctuate by a few cm.The automatic phase detection is able to identify also a wiggling water surface.However, since the GPR waves travel 10 times faster in air than in water and because of the vertical resolution (approx.10 cm) of the GPR, the effect of changing UAS altitude (typically only few cm) on water surface return time is mostly negligible.• A background removal filter was applied to the radargram to increase the signal-to-noise ratio by removing the banded events.The background removal is carried out by subtracting the moving average of the background noise from every trace in the GPR profile.• The signal power of the existing GPR signal was amplified using gains, typically with logarithmic gain function.• The bathymetry layer was extracted by using the layer picking tool provided by GPRSoft with automated phase detection • After the riverbed layer is identified, numeric values, such as trace number, return time and GNSS coordinates, are exported.Furthermore, a file containing NMEA GPGGA messages is exported.This file also contains information about the trace number at which each NMEA sentence (thus each UTC timestamp) was received.• The coordinates in the riverbed layer file are corrected with the PPK position computed in RTKLIB.The PPK-computed positions are matched with the observations of the riverbed layer file by using the relationship between UTC timestamp and trace numbers.• Return time of riverbed is converted into precise depth using the GPR wave speed in water (v water ) computed with Eq. ( 14), in which c is the speed in vacuum and ε water is the real component of relative permittivity of water.The variable ε water was estimated as a function of water temperature, with the relationship shown in Malmberg and Maryott (1956).

Sonar processing
The Deeper Smart Sonar CHIRP+ provides NMEA 0183 messages with depth and position from the in-built GNSS receiver.Since the inbuilt GNSS receiver cannot provide very accurate positions, the UTC timestamps contained in the NMEA 0183 are used to match the PPKcomputed position so that the PPK-computed position can be used to accurately re-estimate sonar position.

Table 2
Site characteristics.table shows water body bottom type and vegetation characteristics, the static electric conductivity of water (σ water ), water temperature (T), coordinates of the markers located on the left and right side of the cross-section (reference is ETRS89-UTM32 for all) and water surface elevation (WSE) in meters above DVR90 geoid.Surveying instrument: "GPR_DB" (drone-borne GPR), "S_100" for sonar 100 kHz, "S_240" for sonar 240 kHz, "S_675" for sonar 675 kHz."RTK_GNSS" (RTK GNSS) and "GPR_WC" (water-coupled GPR) were performed in all cross-sections.

Conversion from depth into bathymetry
Water depth observations, retrieved by sonar or GPR, are subtracted from WSE, retrieved by RTK GNSS, to compute bathymetry, thus elevation of the bed of the water body referenced to the chosen datum (in this research DVR90 geoid).
Lastly, the bathymetric measurements need to be assigned to a position along the cross-section reference line, which is defined by connecting the two reference points on the left and right streambanks.The cross-section line is thus divided into intervals of typically 0.2 m: to each interval we assign the median of bathymetry measurements of each instrument (GPR, sonar or RTK GNSS) that are contained in that interval.Only observations with a maximum perpendicular distance from the observation to the cross-section line of typically 4 m are included in the median.The search radius of 4 m is chosen because the horizontal navigation accuracy of the UAS is not always optimal and an offset of few meters with respect to the cross-section line may occur.

Error statistics
Error statistics, such as Mean Bias Error (MBE), Mean Absolute Error (MAE) and Root Mean Square Error (RMSE), were computed by comparing each RTK GNSS depth measurements to the median of depth measurements (sonar or GPR) contained in the same cross-section discretization interval, if observations were available in that interval.Statistics were computed on depth (not on bathymetry), e.g. a positive MBE for GPR means that the GPR overestimates depth compared to RTK GNSS.

Study areas
The study area consists of two rivers (Guden Å and its tributary Nørreå) and a lake (Furesø), all of which are located in Denmark, as shown in Fig. 3.
Twelve cross-sections in both Guden Å, two in Nørreå, and one lake profile of approx.70 m in Furesø were measured.Table 2 shows the characteristics for the sites, including water temperature and electric conductivity (measured with a commercial electric conductivity meter), bottom type and vegetation density, coordinates of markers located on left and right sides of the profiles and WSE (measured with RTK GNSS).
Furesø presented a sandy bottom, while Guden Å presented a hard silty clay (but with several soft clay patches), and Nørreå had a soft clay with high water content.Vegetation was absent in Furesø.In Guden Å, aquatic vegetation was present in all cross-sections but especially dense in some portions (typically the areas closer to the streambanks) of some cross-sections.The dominant plant species were Batrachium and Glyceria maxima, with the height of the vegetation typically 0.4-0.7 times the water depth.Similarly, in Nørreå, vegetation was sparser than in Guden Å but still present.
Water-coupled GPR was performed in all sites, while sonar and drone-borne GPR were performed in most sites but not all, depending on UAS availability.

Results
In Fig. 4 and Fig. 5 we show the radargrams of four cross-sections (3 in Guden Å and 1 in Furesø), while Fig. 6 compares the bathymetry results of GPR, sonar (three frequencies) and ground-truth RTK GNSS.The radargrams and the bathymetry comparison plots of all the other cross-sections are reported in the Supplementary Material.In the  radargrams, the blue line shows the identified water surface, while the red line highlights the GPR return from the water body bed.The red line is typically drawn by GPRsoft with "automated phase detection" method, while manual correction is performed by the user only in the small portions where automated detection fails (typically close to the streambanks).The radargrams show two significant differences when the antenna is drone-borne (left panels) compared to when it is watercoupled (right panel): i) the electric field-strength of the return from the bed of the water body is weaker and ii) significantly higher noise level in the early portion of the radargram (visible as the white/black bands with strong electric field-strength that propagate from the water surface until approx.70-100 ns).The noise affecting those earlier times is presumably caused by surface cluttering and multipath; indeed, this noise level was not detected when the handheld GPR was elevated (Fig. 1): the UAS carbon fibre structure is suspected to create multiple reflections of the GPR wave.Furthermore, compared to the handheld GPR case, we may assume that the roughness on the water surface induced by the UAS propellers can affect the water surface reflection loss, change wave coherence and increase surface clutter (Arcone and Calkins, 1990;Comite et al., 2017;Rappaport et al., 2003).The case of Furesø, which presents minor bed elevation changes along the crosssection (as visible in Fig. 6), is representative of the surface clutter issue: the noise completely covers the bottom return until the lake is at least 0.8 m deep when the GPR is drone-borne, whilst when the GPR is water-coupled also a depth as shallow as 0.4 m can be detected.
However, in Guden Å cross-sections, also water-coupled GPR is unable to measure close to the stream banks (a.1-2 m horizontal distance from stream edges are unmeasured) for two main reasons: i) the GPR onboard the inflatable boat cannot be dragged up to the exact land-water edge ii) the large GPR footprint is not ideal for measuring the streambanks with steep slopes, causing larger errors in the measured depth close to the streambanks (e.g.due to strong reflections from the dry part) compared to deeper water.As for the Furesø case, in Guden Å the droneborne GPR has also a minimum depth that is site-dependant but is typically in the order of 1 m (e.g.Guden_13859, Guden_13770 and Guden_9158 in Fig. 6) and thus higher than for the water-coupled case.The minimum detectable depth can vary depending on factors, such as i) slope of the water body bed and ii) distance from the streambanks.In terms of riverbed slope, typically higher slope improves the visibility of the returns from the riverbed in the radargram, especially after background noise removal, which typically removes flat horizontal bands such as the lines caused by surface clutter.In terms of distance from streambanks, portions of the water body that are the furthest away from the streambanks are easier to measure because the dry portion of the streambank creates returns that can hide the submerged portion.This is especially critical in light of the fact that low-frequency antennas, such as the Gekko-80, are unshielded and thus can detect radiation reflected from any direction.
The drone-borne GPR profile in Guden_9158 shows a deeper area in the centre of the river (between 17 and 21 m from left bank) that is not visible in water-coupled GPR and RTK GNSS.Among all case studies, this is the only example of significant discrepancy between watercoupled and drone-borne GPR measurements; however, because the significant strength of the signal of both the water-coupled and the drone-borne GPR radargram in that deep portion, it is assumed to be caused by the inaccuracy of UAS navigation: the UAS was flown approx.3-4 m horizontal distance from the cross-section line and thus the error in drone-borne GPR may be associated with local riverbed heterogeneity along the river course.
In all the drone-borne GPR surveys, the GPR antenna was flown at a maximum height of approx.0.5-0.6 m, because higher elevations showed higher noise level, weaker bottom reflection and stronger reflection from the dry portions of the streambanks, making identification of the water body bed more complex.
Fig. 4 and Fig. 5 show only one radargram for drone-borne and water-coupled GPR profile for each cross-section, however multiple profiles were acquired in each cross-sections: radargrams acquired by GPR at the same cross-section were nearly identical, thus GPR measurements are shown to be highly repeatable.
From Fig. 6, it is clear that the sonar measurements at 100 kHz in Furesø are faulty: this is most likely caused by the sediment type and especially by too shallow depth for the 100 kHz frequency (despite a minimum depth of 0.8 m for the 100 kHz, the manufacturer suggests to use this frequency only for depths above 1.5 m).On the other hand, the survey in Furesø highlights that the minimum depth for 240 kHz is approx.0.6 m.
Compared to sonar measurements, the main advantage of GPR is that the GPR can sense through dense aquatic vegetation.This is especially evident for the centre of the cross-sections in Guden_9021 and for the portions close to the streambanks in Guden_13770 (Fig. 6 (b)), Guden_14963 and Guden_8370, in the areas which had high vegetation density according to the visual assessment reported in Table 2.In those cross-sections no sonar frequency (including the 100 kHz) was able to precisely detect the river bed in densely vegetated areas.Guden_8370, shown in Fig. 7, demonstrates that sonar bathymetric measurements show an elevation higher than RTK GNSS, caused by the sonar wave reflection from the vegetation canopy.Indeed, Guden_8370 was highly vegetated on the right side, as it was also shown in Fig. 2, in which the F. Bandini et al. inflatable boat is standing on vegetation.Conversely, drone-borne and water-coupled GPR were not affected by vegetation and were able to penetrate through the submerged canopy.
The statistics comparing all the depth measurements are shown in Table 3; each statistic is computed by comparing the depth measurements of either sonar or GPR with RTK GNSS.Table 3 shows that watercoupled GPR shows slightly better MAE compared to drone-borne GPR in all cross-sections (apart from Guden_13859 and Guden_8828), however the main advantage of water-coupled GPR is that it can measure more shallow water and closer to the streambanks.Average statistics are significantly better for GPR (water-coupled and drone-borne) than for any sonar frequency.
Among the sonar frequencies, the 100 kHz shows the worst performance, as expected from the worst resolution and because it typically penetrates below soft bottoms, whilst the 675 kHz shows the best performance, apart from Guden_13770 (240 kHz is slightly better), Guden_14963 and Nørreå_43716.However, in Guden_14963 and Nørreå_43716 this is an artifact caused by the fact that the 675 kHz was unable (because of connectivity issues) to acquire fully continuous observations in the centre of the cross-section: Thus the statistics of the diffferent sonar frequencies are not comparable for those two crosssections: indeed, the sample size for the statistics computation is smaller for 675 kHz compared to the other two frequencies.
The last row of Table 3 reports the average statistics among all sites; however, direct comparison of the averaged statistics of different instruments should be performed with caution as sonar or drone-borne GPR measurements are missing in some sites.The averaged statistics show that a MAE of approx.8 cm, a MBE of approx. 2 cm and a RMSE of approx.10 cm are achievable with drone-borne GPR, whilst similar statistics are achieved with water-coupled GPR (but on a larger number of sites).
In some cross-sections (e.g.Guden_8917, Guden_9021 and Guden_9158 for water-coupled GPR), RMSE is significantly worse than MAE: this is because RMSE is more sensitive to outliers, specifically large errors occurring in the steep streambanks close to the water's edge (in locations where GPR is reaching the minimum depth capability).
The cross-section in Nørreå_43716 (measured with water-coupled GPR and sonar) shows the worst statistics in terms of MAE and MBE for GPR, which may be caused by overestimation of depth in this crosssection that presented soft bottom.
Among the sonar frequencies, 675 kHz and 240 kHz show negative bias, which translates into underestimation of water depth (vegetation effect), whilst 100 kHz shows positive bias, which may indicate penetration below the bed of the water body.

Discussion
We demonstrated that drone-borne GPR measurements are feasible in waterbodies with a conductivity up to 340 μS/cm and a depth up to 2.5 m.GPR shows slightly worse statistics (depth overestimation) in a cross-section with soft bottom when compared to RTK GNSS.GPR was not tested in waterbodies with higher conductivity or higher water depths, thus the possibility of deploying drone-borne GPR measurements in those conditions remains unknown, especially because the propagation loss increases with electric conductivity and with depth.
Our research shows that GPR measurements give significantly better results than sonar in measuring bathymetry when compared to groundtruth.Ground-truth was taken with RTK GNSS, which has a typical nominal vertical accuracy of approx.2-5 cm; furthermore, when measuring soft riverbeds, the accuracy degrades because it is difficult to keep the bottom plate of the RTK GNSS pole exactly at the bed level instead of penetrating a few centimetres into the sediment.Furthermore, stones and local heterogeneities can complicate direct comparison between RTK GNSS and the footprint of GPR (or sonar).Additionally, Nobes et al. (2018) found that also suspended sediments can have an influence on the accuracy of GPR bathymetric measurements.For these reasons, the achieved MAE of approx.8 cm for GPR (both drone-borne and water-coupled) when compared to RTK GNSS is considered satisfactory.However, drone-borne GPR requires a minimum water depth of ca. 1 m and shows worst results close to the streambanks, thus is not applicable in narrow and very shallow rivers.Furthermore, the low flight altitude requirement for drone-borne GPR can also significantly complicate navigation in narrow rivers with riparian vegetation.
RTK GNSS bathymetric measurements typically have a degree of subjectivity when the operator chooses the exact point where to measure (e.g.avoiding to measure on stones or taking measurements on stones), thus different surveys can result in slightly different results (unless measurements are taken at the exact same cross-section location).On the other hand, GPR surveys have always shown a high degree of repeatability, thus multiple profiles retrieved along the same crosssections give nearly identical radargrams.However, a small degree of subjectivity is left to the user in choosing i) the return from the water surface and ii) the first point on the return from the bottom, especially in portions of the radargram where the "automated phase detection" method fails.Thus, different experienced users could extract slightly different depth observations (typically in the range of 5-10 cm) from the same radargram, especially in the areas close to the streambanks.Regarding sonar deployment, future research could investigate the usage of a higher frequency sonar (above 1 MHz) to have a better accuracy at shallow depths.Furthermore, the deployed sonar provides water depth measurements in real-time from automatic full waveform interpretation; however, if the user had the possibility to extract the full waveform as tabulated data, enhanced removal of vegetation and better interpretation of the bottom return from the full waveform data could improve results in post-processing.Finally, integration of drone-borne GPR with techniques such as drone-borne topographic LiDAR (or photogrammetry) could be performed for estimating the elevation of both submerged (GPR) and dry (topographic LiDAR) portions of the cross-sections, while interpolation could be performed in the areas between dry and submerged parts.

Conclusions
This study assessed the performance of a water-coupled and droneborne GPR for mapping bathymetry in inland water bodies.Significant signal losses (especially due to water surface reflection and spreading through air) affect elevated GPR antennas compared to water-coupled GPR: these losses were computed and subsequently demonstrated with an experiment in which a handheld antenna (elevated at approx. 1 m above the water surface) was compared to a water-coupled GPR.
Field measurements of bathymetry were taken in three different sites (total surveyed cross-sections: 15 of which 1 in a lake) with drone-borne GPR, water-coupled GPR data, sonar and ground-truth RTK GNSS.
The following results were obtained: • Both water-coupled and drone-borne GPR measurements (with GPR antenna flown at approx.0.5 m above the water surface) could detect the bed of the waterbodies.• Two factors are required to effectively map bathymetry with droneborne GPR; (i) a maximum GPR elevation above the water surface of approx.0.5 m and (ii) a minimum water depth in the range 0.8-1.1 m depending on the inclination of the riverbed and distance to streambanks.• Compared to drone-borne GPR, water-coupled GPR can measure shallower water depth (as low as 0.3-0.4m), as the noise caused by surface clutter on the early portions of the radargram does not occur when the GPR antenna is water-coupled.• Both water-coupled and drone-borne GPR provided accurate results with an average MAE of approx.8 cm (thus similar to the vertical resolution of the GPR at this frequency) compared to in situ RTK GNSS.No significant bias was identified in GPR measurements compared to RTK GNSS.• The applicability of the GPR bathymetry monitoring solution was not limited by the conditions at the chosen survey sites in relation to aquatic vegetation, water electric conductivity (tested up to 340 µS/ cm), or maximum water depth (tested up to 2.5 m).GPR was applicable also in waterbodies with soft clay bed, but resulted in slightly worse statistics when compared to ground-truth.the different losses are shown in Table A1.Some of these variables are sensitive to temperature.For example, ε water varies significantly with water temperature from ≈80 at 20 • to ≈87.7 at 0 • (Andryieuski et al., 2015;Malmberg and Maryott, 1956); similarly, a 2 % increase of σ water per 1 • C increase of water temperature typically occurs (Hayashi, 2004).
The variables ε water and ε '' water also depend on the angular frequency of the GPR wave (Andryieuski et al., 2015); however, for frequencies below 1 GHz (as for the GPRs used for bathymetry), the frequency effect on ε water is not significant.
In the Supplementary material, an example of calculations of the different loss coefficients is documented as a Jupyter Notebook.For example, using the variable values from Table A1, the losses of Table A2 were computed.

Table A2
Loss coefficients for each study area, for water-coupled or drone-borne GPR antenna.All units in decibels.

Fig. 1 .
Fig. 1.(a) Different loss components for GPR propagation through air, through water and reflection from the riverbed; (b) radargram showing the returns from the lakebed (approx. 1 m water depth) acquired by a GPR antenna (80 MHz) water-coupled (trace number between 0 and 180 and after 780) and slowly elevated above the water surface (e.g. at trace 430 antenna was elevated approx. 1 m above the water surface).(c) comparison of the lakebed return electric field-strength detected from water-coupled (trace 20) and elevated (trace 430) GPR antenna: the electric field-strength of the lakebed return is approx.5290 bits for water-coupled, whilst only 1190 bits for the elevated antenna.

Fig. 4 .
Fig. 4. Furesø radargram is shown in (a) water-coupled and (b) drone-borne.Guden_13770 is shown in (c) water-coupled and (d) drone-borne.Color lines indicated interpreted water surface and bottom returns.

F
.Bandini et al.

Fig. 6 .Fig. 7 .
Fig. 6.Plots comparing the bathymetric measurements taken with GNSS, sonar (three different frequencies) and GPR.WSE, as measured by the RTK GNSS, is plotted as a flat blue line.Sites: (a) is Furesø lake (cross section is only a small portion of the ca. 3 km lake extension) (b) Guden_13770 (c) Guden_13859 (d) Guden_9158.All elevations are relative to DVR90, which is the standard vertical reference system (geoid) used in Denmark.(For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Table 3
Statistics comparing the depth measurements retrieved by the different instruments, benchmarked against RTK GNSS.MAE is Mean Absolute Error, MBE is Mean Bias Error, RMSE is Root Mean Square Error.Surveying technique: "GPR_WC" (water-coupled GPR), "GPR_DB" (drone-borne GPR), "S_100" for sonar 100 kHz, "S_240" for sonar 240 kHz, "S_675" for sonar 675 kHz.The acronym n.a.stands for non-available.Last row reports the averaged statistics among all sites.The lowest MAE and RMSE among the different survey techniques was highlighted with bold font.
• The sonar provided unreliable (biased) results in vegetated water bodies, but performed well in non-vegetated water bodies when operated at 240 kHz and 675 kHz.(-)