The laser scanner comprises a frequency-doubled IR laser, emitting pulses with about 1.5 ns pulse duration at a wavelength of 532 nm and at a PRR of 50–200 kHz. At the receiver side, the incoming optical echo signals are converted into an electrical signal, they are amplified and digitized at a digitization rate of close to 2G samples/s. The laser beam divergence can be selected between 1 mrad and 6 mrad in order to be able to maintain a constant energy density on the ground for different flying altitudes and therefore balancing eye-safe operation with spatial resolution. The receiver's iFOV (instantaneous Field of View) can be chosen between 3 mrad and 18 mrad. For topographic measurements and very clear or shallow water, a smaller setting is suitable while for turbid water it is better to increase the receiver's iFOV in order to collect a larger amount of light scattered by the water body.

The VQ-840-G employs a Palmer scanner generating an elliptical scan pattern on the ground. The scan range is ±20◦ across and ±14◦ along the intended flight direction and consequently, the laser beam hits the water surface at an incidence angle with low variation. The scan speed can be set between 10–100 lines/s to generate an even point distribution in the center of the swath. Towards the edge of the swath, where the consecutive lines overlap, an extremely high point density with strongly overlapping footprints is produced.

The onboard distance measurement is based on time-of-flight measurement through online waveform processing of the digitized echo signal [62]. A real-time detection algorithm identifies potential targets within the stream of the digitized waveform and feeds the corresponding sampling values to the signal processing unit, which is capable of performing system response fitting (SRF) in real-time at a rate of up to 2.5M targets per second. These targets are represented by the basic attributes of range, amplitude, reflectance, and pulse shape deviation and are saved to the storage device, which can be the internal SSD a removable CFast© card, or an external data recorder via an optical data transmission cable.

Besides being fed to target detection and online waveform processing, the digitized echo waveforms can also be stored on disc for subsequent off-line full waveform analysis. For every laser shot, sample blocks with a length of up to 75 m are stored unconditionally, i.e., without employing any target detection. This opens up a range of possibilities including pre-detection averaging (waveform stacking), variation of the target detection parameters or algorithm, and employing different sorts of full-waveform analysis. For the work presented here, offline SRF with modified target detection parameters was employed for optimized results. The depth performance of the instrument has been demonstrated to be in the range of 2 Secchi depths for single measurements, i.e., without waveform stacking.

Instead of or in addition to the internal camera, a high-resolution camera can be externally attached and triggered by the instrument. This option was chosen for the experiments in this study using a nadir looking Sony Alpha 7RM3. Figures 1 and 2 show conceptual drawings of the sensor and the general data acquisition concept, respectively. respectively. Table 1 summarizes the main sensor parameters.

**Figure 1.** Computer Aided Design rendering of the VQ-840-G laser scanning system in basic configuration. The electrical and data interfaces are shown on the left hand side while the apertures for the green laser scanner and the camera at the bottom of the unit are shown on the right hand side.

**Table 1.** Sensor key specifications.

**Figure 2.** Visualisation of the VQ-840-G mounted on a drone in operation. The scan pattern on the ground is indicated by the green ellipse. The inserts show typical pulse shapes of the outgoing laser pulse and the echo return when measuring into the water.

#### **4. Materials And Methods**

In this Section, the study area is introduced in Section 4.1, the existing datasets are presented in Section 4.2, and the employed data processing methods are detailed in Section 4.3.

#### *4.1. Study Area*

The study area Neubacher Au (N 48◦12'50", E 15◦22'30"; WGS 84) is located in the eastern part of Austria at the tail-water of the pre-Alpine Pielach River, a right-side tributary of the Danube (cf. Figure 3). The study area is part of the conservation area Niederösterreichische Alpenvorlandflüsse (Area code: AT1219000) of the European Union's Natura2000 program. The river is classified as riffle-pool type [63] with a maximum depth of approximately 3 m and exhibits a pluvio-nival regime with expected discharge peaks during snowmelt in winter/spring and torrential rain in summer. Within the investigated reach, the mean annual discharge is about 7 m3/s, the bed-load sediment is dominated by coarse gravel (2–6.3 cm), and the average gradient is about 0.4% [64]. The entire catchment area measures 590 km2 and the mean channel width is approximately 20 m. Although the longitudinal continuum of the Pielach River is disrupted by weirs built for hydropower use and engineering measures, the river has retained some of its natural self-forming morphological characteristics like periodically inundated sidearms, dynamic gravel bars, large woody debris, small oxbows, etc., within the study area [65].

**Figure 3.** Study area Neubacher Au; (**a**) overview map of Austria with study area marked as a red circle; (**b**) orthophoto of the study area overlaid with flight trajectories (circles) and detailed analysis units (rectangles).

The floodplain to the north of the river features around a dozen groundwater supplied ponds (cf. right side of Figure 3). The bottom of the ponds is covered with a layer of mud, which is constantly stirred up by the fish population (e.g., *Cyprinidae*). This causes turbidity resulting in a Secchi depth of generally less than 2 m. With an overall depth of around 5–6 m, the ponds are ideal water bodies for testing the depth performance of bathymetric LiDAR sensors with a nominal depth performance of 1.5–2 SD.

Since 2013, the study area has been repeatedly captured with various topo-bathymetric laser scanners mounted on manned aircraft for monitoring river morphodynamics and instream habitats [4] as well as with UAV-borne topographic sensors for detailed analysis of the riparian area [66,67]. At the same site, first performance evaluations of a UAV-borne bathymetric laser range finder were carried out [61]. All these articles, as well as the cited literature therein, provide further details about the study area.

#### *4.2. Datasets*

In order to assess the accuracy potential and depth performance of the novel topo-bathymetric full-waveform laser scanner, two field campaigns were conducted on 28 August 2019, and 4 September 2019. The two dates are referred to as day 1 and day 2 in the following. On both days, the topo-bathymetric UAV-scanner was operated on a RiCOPTER-M octocopter platform and collected data of the ponds (flight lines marked with red dots in Figure 3) and of an approximately 650 m section of the Pielach River (orange dots).

For comparison and evaluation purposes, two additional LiDAR datasets were captured with the topographic UAV laser scanner VUX1-UAV mounted on a RiCOPTER UAV platform on day 2 and with the VQ-880-G sensor installed on a Diamond DA42 aircraft on 3 September 2019. The respective flight lines are plotted in Figure 3 (VUX1-UAV=blue dots, VQ-880-G=black dots). Both datasets served for cross-validation of the topo-bathymetric UAV-LiDAR data.

Both UAV-based sensor systems operate RGB cameras along with the laser scanners. While the VQ-840-G carries a Sony α7R III (ILCE-7RM3) full-frame mirror-less interchangeable lens camera for capturing 42 MPix RGB nadir images, the VUX1-UAV is equipped with two oblique looking 24 MPix Sony α6000 (ILCE-6000) APS-C format cameras matching the 230◦-FoV of the VUX1 sensor. The RGB images of the latter served as basis for the orthophoto shown in Figure 3. The α7R images, in turn, were used for evaluating the visual depth of the different water bodies during LiDAR data acquisition. Due to technical problems, α7R images were only available for the pond flight lines on day 1 but not for the river flight lines.

The turbidity of the Pielach river depends on the flow velocity and discharge within the entire 590 km2 catchment area. Figure 4a shows the annual mean daily discharge together with the lines for mean discharge (MQ = 7.09 m3s−1) and the one-year flood peak (HQ1 = 106.5 m3s−1), respectively. Discharges beyond the MQ-level only occurred during winter and spring, while summer 2019 was extremely dry with discharges around 2.5 m3s−1 throughout the entire summer. In the days before the data acquisitions, typical summer thunderstorms entailed contamination with suspended debris (cf. Figure 4b) and, thus, a high turbidity level of the river water. While this is generally sub-optimal for laser bathymetry, these conditions also provided a chance to study the depth performance of the scanner and the accuracy implications within the available small range of water depths (Pielach river: ≈3 m, ponds: ≈5–6 m).

**Figure 4.** Discharge curves for study area Neubacher AU; (**a**) mean daily discharge for 2019 (blue) together with nominal mean discharge (MQ, red) and one-year flood peak (HQ1, green); (**b**) minimum, mean, and maximum daily discharge within the data acquisition period (flight dates marked with black ticks).

Table 2 summarizes the in-situ SD and turbidimeter measurements. Secchi depths were measured with a 30 cm black-and-white disk and for measurement of the concentration of suspended particles, a Hach 2100Q portable turbidimeter was employed. From both the SD- and turbidimeter values documented in Table 2, it can be stated that (i) pond FP2 tends to be more turbid than FP1; and (ii) turbidity was generally higher for day 2 compared to day 1. This also relates to the higher discharge peak of 15.7 m3s−1 caused by the thunderstorm event on 29 August before the second data acquisition day compared to the lower peak of 10.5 m3s−1 resulting from the rain event on 24 August before day 1 (cf. Figure 4b).

**Table 2.** Secchi depth and turbidity measurements.

During data acquisition, different scanner settings were tested w.r.t. their implications on depth performance, accuracy potential, and spatial resolution. As stated in Section 3, the bathymetric UAV-laser scanner features both a user-definable laser beam divergence and receiver iFoV. The prior is responsible for eye-safe operation and was set to 2 mrad for all tests, while the PRR, flying altitude, and receiver iFoV were varied. Table 3 summarizes all parameter combinations. In total, the entire test procedure entailed six flights on day 1 and three flights on day 2. The overall duration of each flight amounted to 13 minutes including data acquisition, take-off, initialization of the GNSS/IMU navigation system, and landing. The highest possible PRR of 200 kHz could not be used due to the high power consumption of this scan mode.

**Table 3.** Flight mission parameters.

For validation of the scan data, several reference targets were positioned in the pond and river. To evaluate the effective penetration depth of the bathymetric UAV-scanner, four white-colored metal plates of 1 × 1 m2 size were lowered into pond FP1 on day 1. Each plate was floating in an individual constant water depth, with depths varying from 1.1 m to 3.1 m, corresponding to 0.8–2.2 SD. The floating of the plates was realized by attaching ropes to the four corners of each plate and mounting styrofoam buoyancy bodies to the other end of the ropes. Figure 5 shows the four plates in one of the α7R images during data acquisition on day 1.

**Figure 5.** The 1 × 1 m2 plates (P1–P4) floating in pond FP1 for evaluation of depth performance; camera: α7R; acquisition: day 1; depth of plate from left to right: 1.1 m, 1.9 m, 2.6 m, 3.1 m.

To assess the bathymetric accuracy, an array of three parallel profiles with 6–7 checkerboard targets (i.e., black-and-white disks of 20 cm or 30 cm diameter, respectively) were anchored in the gravel riverbed and measured with a Leica TS16 total station on day 1. The approximate position of the reference target field is marked as R1 in Figure 3. The 19 reference targets exhibited a minimum, mean, and maximum depth of 0.33 m, 1.44 m, and 2.15 m, respectively. Snorkeling was employed for fixing the checkerboard targets in the riverbed as well as positioning the measurement pole in the target's center hole. This was especially necessary for all targets featuring a depth of more than 1 m. As a guide for both installation and measurement, a rope was fixed on both sides of the riverbank. The photographs in Figure 6c–e illustrate the setup of the reference target field.

Furthermore, Figure 6a shows one of the black-and-white checkerboard targets used as ground control points (GCP) for absolute orientation of the image block as well as the laser scans. In total, 25 of theses quadratic targets (edge length: 30 cm) were distributed within the entire study area and measured with a Spectra Precision SP80 GNSS receiver in real-time kinematic (RTK) mode based on correction data obtained from the local GNSS service provider EPOSA.  While GCPs are the first choice for absolute orientation of an image block, laser scanning requires planar control patches. Two A-shaped saddle roofs composed of wooden formwork panels, one on each side of the river in close vicinity to the reference target field, served as the basis for the absolute orientation of the laser flight block (cf. Figure 6d,e). The data captured on day 2 were only used for orthophoto generation and cross-validation on dry land due to the very high turbidity on this day. The topo-bathymetric LiDAR data of day 1, in turn, serve as a basis of all subsequently described data processing steps and the presented results.

**Figure 6.** Photo documentation of data acquisition on 27 August (day 1); (**a**) Pielach river with quadratic black-and-white checkerboard ground control point target in the foreground; (**b**) topo-bathymetric UAV laser scanner RIEGL VQ-840-G during data capture; (**c**) UAV-image of field of submerged circular black-and-white checkerboard targets aligned with rope fixed between both riverbanks; (**d**) terrestrial photo of reference target field, most of the targets are not visible with the naked eye due to the limited Secchi depth of 1.1 m, yellow wooden plates in fore- and background: saddle roof control patches for absolute orientation of laser scans; (**e**) total station Leica TS16 in front of reference field and control patches.

#### *4.3. Data Processing Methods*

In the following, the general workflow including preprocessing, refraction correction of the raw LiDAR point clouds, assessment of accuracy and spatial resolution, and evaluation of depth performance is summarized. For most of the processing steps, the scientific laser scanning software OPALS [68] was employed.

#### 4.3.1. Preprocessing

Data preprocessing comprised the following steps:

#### 4.3.2. Evaluation

	- *Measurement noise:* Assessment of the sensor's noise level was carried out at smooth surfaces (e.g., asphalt) within the individual strips by estimating a best-fitting plane for each point of an appropriate surface patch via least-squares adjustment based on the k nearest neighbor points (k = 12) within a maximum search radius of 50 cm. For surfaces with known planar geometry, the residuals (i.e., point-to-plane normal distances) provide an estimate of the sensor's ranging precision (cf. Figure 7a).

#### **5. Results**

In this section, we present the quality evaluation results separated into precision (Section 5.1), accuracy (Section 5.2), depth performance (Section 5.3), and spatial resolution (Section 5.4).

#### *5.1. Precision*

#### 5.1.1. Measurement Noise

Figure 8 shows five selected surface patches (asphalt: A1, A2, gravel: G1, G2, and meadow: M1). For each patch we provide the orthophoto and DEM (color-coded height map superimposed with 1 cm contour lines) in the first two rows. The last row exhibits a plot of the local interpolation error (σ0). In case of a sufficiently planar surface, σ0 quantifies the measurement noise (reproducibility). This is especially the case for A1. All other patches exhibit small scale topographic variations (gravel, G1 and G1) or undulations due to grass vegetation (A2, M1).

**Figure 8.** Visualization of measurement noise for selected 5 × 5 m2 surface patches; A1: asphalt, A2: asphalt and grass, G1/G2: gravel bank, M1: meadow; first row: digital orthophoto, second row: color coded height map superimposed with 1 cm contour lines, third row: color coded map of local surface roughness (i.e., σ0 [m]).

Table 4 summarizes the quantitative results. For the asphalt dominated patches Table 4 also documents the noise level variations depending on the mission parameters (flying altitude, iFoV). For A1, σ0 is below 2 mm for a narrow iFoV of 6 mrad and slightly higher (3.5 mm) for the wider iFoV of 18 mrad. For the mixed asphalt/grass patch A2, increasing the flying altitude from 50 m to 70 m does not effect the precision (mean: ≈4 mm, median: ≈3 mm). In both cases, the higher mean value reflects the higher variations within the grass vegetation, while the more robust median can be seen as representative for the dominant asphalt area. For all other patches, the micro-relief leads to increased σ0 values, which are in the range of 1 cm for the gravel patches (G1, G2) and around 5 mm for the meadow patch M1. All reported values are below 15 mm, the precision value reported in the sensor's spec sheet [9].

**Table 4.** Evaluation of measurement noise reported as local σ0 (mean, median, and interquartile range (IQR)).

#### 5.1.2. Strip Fitting Precision

Figure 9 shows a color-coded map of strip height differences in smooth surface areas on top of a DSM map of the entire UAV-LiDAR flight block (color-coded elevation map superimposed with shaded relief map). The dominant white color tones of the strip differences indicate deviations of less than 2 cm and confirm the good fitting precision of the flight strips. The histogram of the strip height differences plotted in the lower-left corner of Figure 9 exhibits a symmetric shape. The distribution is unbiased (mean: 1 mm, median: 0 mm) with a dispersion of 9 mm (σ*MAD*). Occasional strip differences in vegetated areas, which have passed the automatic smoothness filter, are responsible for the comparably high standard deviation of 10 cm. 96% of the strip height differences are within a range of −24 mm to +25 mm, which further underlines the fidelity of the entire sensor system including navigation device, laser ranging unit, and scanner.

**Figure 9.** Color coded mosaic of strip height differences, background: LiDAR DSM (color coded elevation map superimposed with shaded relief map), lower left: histogram of height differences in smooth surface areas, lower right: color bar for strip differences.

#### *5.2. Accuracy Assessment*

The absolute accuracy, measured as signed point-to-plane distances between terrestrially measured reference points and LiDAR-derived planes are plotted in Figure 9 (labeled red dots) and Figure 10 for land and underwater, respectively.

On the land side, the deviations are mainly below 3 cm, with only two checkpoints at the western end of the flight block exhibiting deviations of ≈5 cm. The poorer absolute fitting accuracy in this part of the block is due to the unavailability of appropriate laser control patches in this area of the block, where absolute orientation was of minor concern for the conducted evaluation. Control patches (wooden saddle roof formwork constructions) were located in the vicinity of the river section featuring the submerged reference targets (cf. Figure 6c–e). In this area utmost accuracy of the absolute block orientation was crucial, and the deviations in this region are 2 cm or better. Table 5 summarize the deviations on dry land for selected reference points on the saddle roof shaped control patches (CP, cf. blue dots in Figure 9) and on the water-land-boundary (WLB, cf. black dots in Figure 9). The deviations are generally below 1 cm (cf. points 1001–1020, 2001–2010, 108–110) and only the points on the ridge of the control patches (1022–1024, 2010–2012) exhibit larger, systematically negative deviations of more than 3 cm. This will be further discussed in Sections 5.4 and 6.

**Table 5.** Accuracy assessment for reference targets on land; selected point-to-plane distances [m] between reference points and LiDAR derived planes.

Figure 10 shows the results of the bathymetric accuracy assessment for all 20 submerged reference targets. The upper left plot shows the locations of the reference targets, arranged in three cross-sections, together with the point labels of the middle section, for which the deviations are presented numerically in Table 6 together with the respective water depth. For refraction correction of the raw laser echoes, we strictly used a refractive index of *n* = 1.33 for beam deflection (i.e., angular component). For run-time correction, we varied the coefficient from *n* = 1.33 to *n* = 1.36 according to literature [75,76] and also beyond that till *n* = 1.42. Table 6 shows that the use of standard refraction coefficients (1.33, 1.36) leads to a systematic, depth-dependent, negative bias indicated by the red color tones. The respective columns in Table 6 also clearly show the depth-dependent increase of the deviations reaching 7.8 cm for the deepest point 206 (water depth = 2.1 m) for refractive index of 1.36. Due to this unexpected systematic effect, we first double-checked the accuracy of the water surface, but did not find any evidence that the modeled water surface model is responsible for the bias (cf. deviations water-land-boundary points in Table 5 and the cross-sectional plot in Figure 11). In addition, we performed refraction correction with different software tools (RiPROCESS, OPALS) leading to consistent results. We, therefore, performed tests with refraction coefficients beyond *n* = 1.36 and found that the bias entirely disappears for *n* = 1.42. This can be seen from the white color tones of the respective plots in Figure 10 as well as from the underlined values in Table 6.

**Figure 10.** Color coded point-to-plane distances [m] between check points and LiDAR derived planes, run-time correction based on refraction coefficients ranging from *n* = 1.33 to *n* = 1.42.

**Table 6.** Accuracy assessment for underwater reference targets; selected point-to-plane distances [m] between reference points and LiDAR derived planes using different refraction coefficients for run-time correction.

Figure 11 shows a vertical section of the refraction corrected LiDAR point cloud (*n* = 1.36) colored by signal amplitude overlaid with (i) the reference points (magenta) and (ii) the modelled water surface (black). The saddle roof on the right (i.e., western) embankment shows very good agreement with the reference points whereas the submerged reference targets lie systematically above the LiDAR point cloud of the river bottom.

**Figure 11.** Vertical section of refraction corrected 3D point cloud (*n* = 1.36) of strip 105334 colored by signal amplitude; larger magenta dots: checkpoints, black: modelled water surface.

#### *5.3. Depth Performance*

To evaluate the sensor's depth performance, four white metal plates were lowered into pond FP1 at the eastern border of the flight block (cf. Figure 5) in depths of 1.1–3.1 m corresponding to 0.8–2.2 SD. The plates were captured multiple times with different sensor settings (PRR, altitude, iFoV). Figure 12 depicts the 3D point clouds of two selected strips (063116, 093054), both flown 50 m above ground level with a PRR of 50 kHz, but with a different iFoV of 6 mrad and 18 mrad, respectively. Figure 12 reveals that the wider iFoV generally delivers more points on the detectable plates P1-P3 and the plates appear larger than the nominal size of 1 m2. The signal amplitude is higher for the wider 18 mrad iFoV, indicated by the more reddish color tones, as more backscattered laser light reaches the detector.

**Figure 12.** Perspective view of 3D point points of lowered metal plates P1-P4 colored by signal amplitude [dB]. upper row: strip 063116, iFoV: 6 mrad, lower row: strip 093054, iFoV: 18 mrad.

For the same strips, Table 7 provides a quantitative summary. Table 7 documents that plates P1–P3 could be detected in the point clouds, but not P4 at a depth of 3.11 m or 2.2 SD, respectively. This holds for the strip flown with a 6 mrad iFoV as well as for the strip with the wider 18 mrad iFoV. As expected, the signal amplitude level of the 18 mrad iFoV is generally higher for all plates, and at P3 the receiver still measures 14.5 dB. This, however, was not enough to identify P4 located just beyond the 2 SD limit. The relative depth of P3 amounts to 1.94 SD, thus, slightly below the sensor specification of 2 SD. Furthermore, the nominal-actual depth comparison reveals the same systematic depth overestimation as already reported in Section 5.2.

**Table 7.** Depth performance assessment.

Figure 13 further illustrates the impact of the user-definable receiver iFoV in terms of depth penetration at a near shore transect of pond FP1. In both cases, the point clouds were captured from a flying altitude of 50 m but with a different iFOV (strip 063419: brown, 6 mrad; strip 093410: blue, 18 mrad). The larger iFoV leads to a higher depth penetration. In the depicted case of a freshwater pond with an SD of 1.4 m at the time of data acquisition, the maximum depth is about 2.67 m with the 6 mrad iFoV and 3.07 m using the larger 18 mrad iFOV, corresponding to 1.9 SD and 2.2 SD, respectively. The gain in the presented case is 40 cm (0.3 SD) or 15%, respectively. Thus, the manufacturer's specification of 2 SD was verified in the experiment. The wider iFoV also marginally increases the measurement noise, which can be seen from the higher amount of clutter points in Figure 13 as well as from Table 4 (A1, median). This aspect is further discussed in Section 5.4.

**Figure 13.** Vertical section of refraction corrected 3D point clouds of strip 063419 (brown, iFoV: 6 mrad) and 093410 (blue, iFoV: 18 mrad); flying altitude: 50 m; refraction coefficient used for run-time correction: (*n* = 1.36).

#### *5.4. Spatial Resolution*

The spatial resolution of any LiDAR system is limited by the size of footprint diameter, which, in turn, is a function of the beam divergence γ and the measurement range *R*. In this experiment, a constant beam divergence of 2 mrad was employed to ensure eye safety at the standard flying altitude of 50 m. The footprint diameter *D* calculates to *D* = *R*γ, i.e., 10 cm for the standard flight setup in our experiment.

To evaluate the actual footprint diameter, two suitable objects have been chosen and analyzed in detail, namely (i) a linear object and (ii) a sharp angle. Figure 14 shows the 3D point clouds of the study objects in a perspective view (a-c) and a sectional view (d) colorized by signal amplitude. The linear object (rope) does not appear as a confined line of points but rather as a band featuring a width of 30 cm, thus, 3 times the theoretical footprint size.

**Figure 14.** Empirical assessment of spatial resolution; top row: rope spanned across river, bottom row: gable roof shaped control patch; (**a**–**c**) perspective view of 3D point cloud; (**d**) vertical section; points colored by signal amplitude [dB].

At the ridgeline of the gable roof-shaped control patch, the points do not form a sharp edge but appear as a circular cylinder featuring a diameter of 20 cm as the finite laser footprint illuminates both roof sides. This effect is responsible for the larger deviations for all control points along the ridgeline reported in Section 5.2 (cf. Table 5, e.g., point 1022) while all control points, for which the laser footprint is entirely on one side of the wooden panel, exhibit deviations below 1 cm.

It is noted that, as expected, also the receiver iFoV substantially influences the spatial resolution, even for constant beam divergence settings. The plate experiment presented in Section 5.3 shows this clearly. While all plates exhibit an area of 1 m2, increasing the iFoV has a big impact on the apparent

size of the plate as depicted in the point cloud. For plate P1, for instance, the area covered by laser returns reflecting from the plate with the 18 mrad iFoV setting measures 3.5 m2. This is caused by beam spreading due to volume scattering in the water column. With a larger iFOV, more signal from the periphery of the conical laser beam reaches the detector enabling echo detection even if the plate is only partially within the illuminated footprint. However, full and partial hits of the plates can be distinguished via analysis of the amplitude values. Amplitudes are nearly constant if the entire laser footprint impinges on the plate, and rapidly drop outwards with a decreasing percentage of illuminated plate area (cf. bottom left plot of Figure 12).

#### **6. Discussion**

The presented compact topo-bathymetric laser scanner exhibits a weight of 12 kg and a power consumption of 160 W. This makes the instrument suitable for integration on light manned aircraft, gyrocpters, and helicopters as well as on larger UAV platforms, as was the case for our data acquisition. The considerable payload and power consumption, however, limit UAV flight endurance with the selected UAV. The longest net acquisition airtime in our experiment was about 10 min excluding starting, landing, and IMU/GNSS initialization. This might be sufficient endurance for corridor mapping in a visual line of sight (VLOS) context, but poses hard constraints for larger project areas and beyond VLOS operation. Thus, our UAV-borne experiment is regarded as a successful proof-of-concept, but further sensor miniaturization is needed for integration on mainstream UAV platforms like the *DJI Matrice 600 Pro* or comparable platforms. No matter if the sensor is integrated on a light manned or unmanned aerial platform, in both cases the mobilization costs are considerably lower compared to traditional ALB surveys and smaller aircraft feature a higher maneuverability. The sensor, therefore, constitutes a cost-effective alternative, especially for mapping bathymetry of narrow river channels, as a supplement to area-wide topographic laser scanning which is today often available on a country-wide level.

A major benefit of the sensor design is the adjustable pulse repetition rate, laser beam divergence, and receiver iFoV, making eye-safe operation possible for low flying altitudes like the standard 50 m altitude in our experiment. Due to the smaller sensor-to-target distances and the narrow laser cone, the resulting nominal laser footprint diameter in the dm-range is much smaller compared to existing topo-bathymetric sensors operated on manned aircraft featuring footprint sizes of about half a meter. This gain in spatial resolution opens new applications in the context of hydraulic engineering (roughness estimation) for reliable flood simulations, eco-hydraulics (habitat modelling), fluvial geomorphology, and many more.

As stated before, the limiting factor for spatial resolution in the lateral direction of the laser beam is the *laser footprint*, i.e., the illuminated area on the target. As reported in Section 5.4, this is not the same as the area resulting from the vendor specifications. The power distribution in the laser footprint has Gaussian shape and it is common practice to specify the beam diameter by twice the radius where the optical power has decreased by a factor of 1/*e*2. For a beam diameter *w*, the power distribution *I*(*r*) as a function of off-axis distance *r* is given by

$$I(r) = I\_0 \mathfrak{e}^{-\left(\frac{2\sqrt{2}}{w}r\right)^2},\tag{1}$$

with *I*0 the power at the center. The power level at *r* = *w*/2 is 0.14% of the power level at the maximum. Expressing the power level at *r* = *w*/2 in Decibels gives

$$10\log\_{10}\left(\frac{I(0.5w)}{I\_0}\right) \approx -8.7\text{dB}.\tag{2}$$

A signal causing a return strength in the center of the beam that is more than 8.7 dB above the receiver detection threshold will cause a detectable signal in an adjacent measurement even when further apart than the nominal beam diameter radius. The effective size of the footprint, thus, is determined not only by the properties of the laser beam but also by the sensitivity of the receiver circuitry and the iFOV. A small isolated object will be seen in adjacent laser points until it is too weak to be detected. However, the object will appear with diminishing amplitudes towards the edges (cf. Figure 14b).

The finite footprint and high receiver sensitivity also cause smoothing effects if the target's shape is not flat within the illuminated area. This applies to discontinuities in height (cliff, wall) or slope (sharp edge), or in case of high curvature (e.g., small boulders). In any case, this leads to rounding artifacts as exemplified in Figure 14d for a sharp edge. In the context of mapping river channel bathymetry, this limits the ability to precisely reconstruct objects smaller than ≈20 cm. While the instrument's spatial resolution is too low for roughness estimation on grain size level, it is suitable for the detection of blocks, boulders, and other flow resistance relevant objects. As an example, the circular underwater reference targets clearly stood out from the surrounding river bottom although they were firmly anchored in the ground and, thus, only protruded from the ground surface by a few cm.

For smooth surfaces, the sensor showed a measurement noise below 1 cm verified at horizontal asphalt surfaces (cf. surface patch A1 in Figure 8 and Table 4) as well as on tilted surfaces (cf. Table 5, e.g., points 1001–1020). Beyond the local level, this also holds for the strip fitting precision of the entire flight block, which even outperformed comparable data acquisitions with a state-of-the-art topographic UAV laser scanner employed in the same study area in the previous work [66,71]. High relative accuracy is of special importance, e.g., for detailed hydrodynamic-numerical modeling in flat areas where a few cm can determine whether or not certain areas are flooded.

While the absolute accuracy of the flight block, measured as point-to-plane distances between terrestrially surveyed reference targets and the LiDAR point cloud, is generally better than 3 cm on land (cf. red dots in Figure 9), systematic deviations were observed for submerged checkpoints. In addition to geometric calibration of the sensor system and (ranging, scanning, lever arms, boresight angles, trajectory), the overall accuracy of LiDAR-derived bathymetry depends on the accuracy of (i) the modeled water surface and (ii) the refraction coefficients for run-time correction. Data analysis has revealed that laser echoes from the water surface are dense (i.e., practically every laser shot returned an echo from the surface) and the signal amplitude is sufficiently strong to enable reliable range measurement. Figures 11 and 14a,b illustrate that (i) the recorded signal amplitudes are around 20 dB; (ii) the point density is high; and (iii) the vertical spread is remarkably low. It is well known from literature including our own work [1,57,77,78] that the laser signal reflected from water surfaces is a mixture of reflections from the surface and volume backscattering from the water column beneath the surface. Surface data obtained from conventional green-only ALB sensors operated from manned platforms in flying altitudes of around 600 m often show a fuzzy appearance. In contrast, surface data acquired in this experiment are concise and the point clouds show a vertical spread of only a few cm. Independent validation of the modeled water surface with terrestrially surveyed checkpoints at the water-land-boundary (cf. Table 5) also confirmed the accuracy of the water surface model. Thus, we rule out that inaccuracies of the modeled water surface are responsible for the detected depth-dependent bias.

Concerning run-time correction, we first started with standard values documented in the literature. [76] reports a valid range of *n* = 1.329–1.367, but the upper bound is only valid for extreme salinity and temperature conditions. In recent work, [75] reported that group velocity is relevant for run-time correction of pulsed lasers rather than the phase velocity of the underlying laser wavelength (λ = 532 nm). With a respective experiment, they showed that a value of around 1.36 is valid for freshwater conditions. Therefore, this value was also used for run-time correction in this work, but comparison with the reference points still showed a systematic, depth-dependent error. The bias only disappeared when using a hypothetical refraction coefficient of *n* = 1.42. As there is no physical explanation for such a high refractive index, the question for the reason of the systematic deviation still remains. Multi-path effects due to complex forward scattering within the water columns are the most likely explanation for the observed phenomenon. As documented in Table 2 and as can be

seen from the on-site photographs in Figure 6, the turbidity level was high (Secchi depth in the river: 1.1 m) due to a thunderstorm event in the days before data acquisition. This resulted in a high load of suspended sediment. Scattering at the sediment particles leads to a widening of the laser footprint, which entails non-linearity of the ray paths leading to elongated path lengths in water, which might result in the observed underestimation of the river bottom height. Following this line of argument, a repeat survey is planned in the winter season 2020 in clear water conditions to verify and quantify the influence of turbidity on the UAV-LiDAR derived depth measurements.

The depth performance of 2 SD stated by the manufacturer was confirmed in the conducted evaluation. All results relied on online waveform processing [62] and did not make use of sophisticated waveform processing [79,80]. Processing of data acquired with different iFoV settings confirmed an increased depth performance using the widest possible setting of 18 mrad compared to the default setting of 6 mrad. The gain measures 0.3 SD or 15% but comes at the price of a reduced spatial resolution as the receiver's perceivable target area is larger, exemplary illustrated in Figure 12 at plate P1. In a practical context, this clearly shows that a trade-off between high spatial resolution and maximum depth performance needs to be found. A potential best practice procedure may include a repeat survey of the same flight lines first with scanner settings optimized for spatial resolution (high PRR, narrow beam divergence and iFoV) and second with parameters optimized for maximum depth penetration. With the studied instrument, this can be achieved in a single flight mission as individual sensor settings are feasible on a per flight line basis.

The conducted depth performance assessment based on the four lowered metal plates yielded that plate P4 at a depth of 2.2 SD could not be identified in any of the flight strips including the ones with 18 mrad iFoV. This is astonishing, as the natural bottom of the pond could be measured at approximately the same depth (cf. blue points in Figure 13). By extrapolating the signal amplitude decrease from plate P3 (14.5 dB), enough signal should have been available to also detect P4. One of the plausible reasons is an abrupt turbidity increase within the freshwater pond at a water depth of around 3 m. We consider this the most likely explanation as the pond is mainly used for fishing purposes and the dominating species (*Cyprinus carpio*) tends to dig into the mud at the bottom of the pond potentially swirling up mud particles.

Beyond the pure sensor settings, an additional increase of depth performance is expected via sophisticated waveform processing. [81] reported an increased areal coverage when employing a surface-volume-bottom based exponential decomposition approach, in both very shallow and deeper water areas. [5,82] use waveform stacking, i.e., averaging of neighboring waveforms, to enhance ALB depth performance reporting an increase of 30%. The presented sensor is specifically well suited for the application of waveform stacking, also referred to as pre-detection averaging, as waveform recording does not depend on a trigger event (i.e., echo detection) during data acquisition but waveforms are continuously stored within a certain range gate. While this is not possible for typical ALB flying altitudes of 600 m due to storage capacity limits, it is well possible for the shorter ranges in UAV-borne laser bathymetry. By combining the possibilities offered by the sensor (PRR, beam divergence, iFoV) with data processing based on the cited strategies, a maximum depth penetration of about 2.9 SD can be expected, but further experiments are required to verify this value.

Further advantages arise from the integration of active and passive sensors in a single hybrid system. Examples for applications benefiting from the concurrent acquisition of laser and image data are water body detection [83] and discharge estimation [46].

#### **7. Conclusions**

In this paper, we introduced the concept of a novel compact topo-bathymetric laser scanner and presented a comprehensive assessment evaluating measurement noise, strip fitting precision, accuracy compared to terrestrially measured reference data, depth performance, and spatial resolution. With a weight of 12 kg the sensor is suitable for integration on light manned as well as powerful unmanned platforms, reducing the mobilization costs compared to conventional ALB. The scanner

features an adjustable pulse repetition rate (50–200 kHz), scan rate (10–100 lines/s), beam divergence (1–6 mrad), and receiver's iFoV (of 3–18 mrad). With nominal flying altitudes from 50–150 m, this enables flexible flight planning and full user control concerning the properties of the resulting 3D point clouds w.r.t. point density, spatial resolution, achievable depth performance, and aerial coverage.

For assessing the overall quality, a flight mission was designed and carried out on 28 August 2019, at the pre-Alpine Pielach River (Austria) with the sensor mounted on an octocopter UAV platform. After geometric calibration and refraction correction, the resulting 3D point clouds were compared to reference points obtained by terrestrial surveys with RTK GNSS and total station. The evaluations exhibited a local measurement noise at smooth asphalt surfaces of 1–3 mm, a relative strip fitting precision of about 1 cm, and an absolute flight block accuracy of 2–3 cm compared to check points on dry land.

Assessment of the bathymetric accuracy yielded a depth-dependent bias when employing a representative refraction coefficient of *n* = 1.36 for run-time correction in water. The maximum deviation of 7.8 cm at a water depth of 2.1 m only disappeared when using a hypothetical refraction coefficient of *n* = 1.42, which is physically implausible. We hypothesize that the bias results from multi-path effects caused by forward scattering at dissolved sediment particles, as turbidity was high (SD: ≈1.1 m) due to thunderstorm events in the days before data acquisition.

Depth performance evaluation, based on four metal plates lowered into a freshwater pond in different depths, confirmed a maximum water penetration depth of 2 SD for the laser echoes derived with online waveform processing. In addition to that, the entire echo waveform is also stored for off-line waveform analysis. Waveform recording does not depend on a triggering event but the entire waveform is rather captured within a user-definable range gate. This enables waveform stacking by summing up the waveforms of neighboring laser shots, which potentially increases the depth performance to 2.5–2.8 SD. Respective experiments addressing these open questions are currently in preparation.

To sum up, the presented compact topo-bathymetric laser scanner is well suited for mapping river channel bathymetry. The sensor system poses an alternative to conventional ALB for mapping smaller rivers and shallow lakes when mounted on flexible UAV-platforms and also for larger coastal environments when integrated on light manned aircraft. One of the main benefits compared to other UAV-based bathymetric laser sensors is the full adjustability of the sensor parameters enabling the end-user to balance accuracy, depth performance, spatial resolution, and aerial coverage.

**Author Contributions:** G.M. designed the data acquisition (selection of study area, flight mission parameters, reference measurements) and was responsible for data processing (geo-referencing, performance assessment), writing main parts of the text, and supervision the overall manuscript preparation. M.P. provided the text for sensor concept section. R.S. has performed the waveform analysis and other data processing tasks. S.F. and L.N. designed and installed the measurement field for accuracy assessment, and conducted the terrestrial surveying. All authors contributed to the interpretation of the results and to revision of the manuscript. All authors have read and agreed to the published version of the manuscript.

**Funding:** The presented investigation was partially conducted within the project "Bathymetry by Fusion of Airborne Laser Scanning and Multi-Spectral Aerial Imagery" (SO 935/6-2) funded by the German Research Foundation (DFG). The work of Martin Pfennigbauer and Roland Schwarz was funded by the Austrian Research Funding Agency (FFG) within the project "Aerial search & Rescue support and supErvision of inAccessible terrainS" (AREAS) 86702.

**Acknowledgments:** The authors acknowledge TU Wien Bibliothek for financial support through its Open Access Funding Program.

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

#### **References**

1. Guenther, G.C.; Cunningham, A.G.; Laroque, P.E.; Reid, D.J. Meeting the accuracy challenge in airborne lidar bathymetry. In Proceedings of the 20th EARSeL Symposium: Workshop on Lidar Remote Sensing of Land and Sea, Dresden, Germany, 16–17 June 2000.

© 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* **Quantifying Below-Water Fluvial Geomorphic Change: The Implications of Refraction Correction, Water Surface Elevations, and Spatially Variable Error**

#### **Amy S. Woodget 1,\*, James T. Dietrich 2 and Robin T. Wilson 3**

Received: 12 September 2019; Accepted: 15 October 2019; Published: 18 October 2019

**Abstract:** Much of the geomorphic work of rivers occurs underwater. As a result, high resolution quantification of geomorphic change in these submerged areas is important. Currently, to quantify this change, multiple methods are required to get high resolution data for both the exposed and submerged areas. Remote sensing methods are often limited to the exposed areas due to the challenges imposed by the water, and those remote sensing methods for below the water surface require the collection of extensive calibration data in-channel, which is time-consuming, labour-intensive, and sometimes prohibitive in difficult-to-access areas. Within this paper, we pioneer a novel approach for quantifying above- and below-water geomorphic change using Structure-from-Motion photogrammetry and investigate the implications of water surface elevations, refraction correction measures, and the spatial variability of topographic errors. We use two epochs of imagery from a site on the River Teme, Herefordshire, UK, collected using a remotely piloted aircraft system (RPAS) and processed using Structure-from-Motion (SfM) photogrammetry. For the first time, we show that: (1) Quantification of submerged geomorphic change to levels of accuracy commensurate with exposed areas is possible without the need for calibration data or a different method from exposed areas; (2) there is minimal difference in results produced by different refraction correction procedures using predominantly nadir imagery (small angle vs. multi-view), allowing users a choice of software packages/processing complexity; (3) improvements to our estimations of water surface elevations are critical for accurate topographic estimation in submerged areas and can reduce mean elevation error by up to 73%; and (4) we can use machine learning, in the form of multiple linear regressions, and a Gaussian Naïve Bayes classifier, based on the relationship between error and 11 independent variables, to generate a high resolution, spatially continuous model of geomorphic change in submerged areas, constrained by spatially variable error estimates. Our multiple regression model is capable of explaining up to 54% of magnitude and direction of topographic error, with accuracies of less than 0.04 m. With on-going testing and improvements, this machine learning approach has potential for routine application in spatially variable error estimation within the RPAS–SfM workflow.

**Keywords:** fluvial; geomorphology; change detection; remotely piloted aircraft system; refraction correction; structure-from-motion photogrammetry; water surface elevation; topographic error; machine learning

#### **1. Introduction**

High spatial resolution, spatially continuous, accurate, and precise methods for quantifying topographic change in fluvial environments are important for a range of river science and management applications. For instance, such data can provide valuable input to sediment routing and budgeting models (e.g., [1,2]), help predict channel dynamics through the provision of channel texture and roughness data [3–6], assist in the monitoring of patterns of geomorphic adjustment [7,8] to significant external drivers such as climate change [9], and reveal trends in channel capacity, which in turn feed flood risk models and insurance forecasting [10,11]. Such data are also essential in assessments of habitat quality and availability [12,13]. Significant early work on quantifying topography and geomorphic change, some of which has been undertaken in fluvial environments, made use of photogrammetry [14,15], terrestrial and airborne laser altimetry [8,16,17], and extensive GPS- or total station-based topographic surveys [7,18]. Such studies tend to use either multiple approaches for characterising the exposed and submerged parts of the system separately, which are subsequently fused together, or dense networks of irregularly spaced individual point measurements, which are exceptionally time-consuming to collect and must then be interpolated to form an elevation model. The resulting models typically have spatial resolutions no finer than 1 m/pixel. Ideally, a single method of quantifying topographic change in all parts of the river system could be used to enable more rapid and consistent data acquisition for monitoring and management of river habitats. Such a method should have high accuracy and precision, a very high spatial resolution (i.e., at the grain scale, <0.1 m/pixel), and provide spatially continuous data (typically as rasters or point clouds) rather than interpolations between sparse point samples (such as traditional topographic surveys, total station surveys, global navigation satellite system (GNSS) surveys, and cross sections or longitudinal profiles). This ideal method would be in alignment with the popular 'riverscape' paradigm [19], which conceptualises fluvial systems as highly heterogeneous rather than gradually varying and would provide physical habitat data at the grain scale, which is of most relevance to aquatic biota (i.e., <0.1 m spatial resolution over reach-scale extents).

Recently, we have seen a proliferation of studies offering high spatial resolution topographic data over continuous spatial extents in a variety of geomorphological settings, without the need for interpolation between data points (e.g., [4,8,20–26]). Such studies have demonstrated the use of airborne and terrestrial LiDAR and Structure-from-Motion photogrammetry using imagery from remotely piloted aircraft systems (RPAS). In fluvial environments specifically, the acquisition of survey data is complicated in submerged areas where the effects of refraction, turbidity, and water surface roughness represent well-known challenges [27,28]. However, these submerged spaces are where a majority of fluvial 'work' occurs (i.e., erosion and deposition resulting from variations in flow magnitude and shear stress). To date, most fluvial studies aimed at measuring spatially-continuous topography and/or topographic change at high spatial resolutions offer just one of the following:

bathymetric approach on images acquired by RPAS (remotely piloted aircraft system) to obtain data in submerged areas. This approach requires the associated acquisition of bathymetry elevations using a GNSS to calibrate the relationship between the spectral image data and the water depth (from which elevation can be inferred). Similar workflows are employed by [38].

As such, no study has yet demonstrated and rigorously assessed an approach which uses a single method capable of quantifying spatially continuous fluvial geomorphic change at grain scales, both above and below the water surface, without the need to collect extensive calibration data. It is here that this paper aims to contribute.

#### *1.1. Refraction Correction*

Since 2015, a method which uses high resolution RPAS (or UAS/unmanned aerial system) imagery, processed using SfM, has been shown as a promising tool for measuring both above-water (bank and floodplain) and below-water (bathymetric) topography at the grain scale in fluvial environments, when paired with refraction correction post-processing [27]. Above-water use of the RPAS–SfM approach has already provided an important step-change in our ability to quantify exposed topography without the greater expense of a terrestrial laser scanning system (e.g., [34,38–40]), but its use in submerged areas is less well established, due to the additional challenges posed by the presence of water and, in particular, the effects of refraction. Recently, two key RPAS–SfM methods have been proposed by the authors for use in submerged areas [27,28].

In 2015, Woodget et al. [27] published details of a simple approach for refraction correction of bathymetric elevation models using a small angle approximation. This RPAS–SfM approach generated high resolution orthophotos and digital elevation models (DEMs) of submerged areas. These datasets were analysed in a GIS environment to extract the position and elevation of the water's edge and create an elevation model for the position of the water surface. Snell's Law was then applied to this surface to increase the apparent water depths by a multiplier of 1.34 (i.e., the refractive index of clear water) to account for the effects of refraction (which cause the channel bed to appear closer to the camera than it really is). This workflow provides a relatively simple method for use in shallow (<1 m), non-turbid settings and is easily implemented by other researchers and practitioners (e.g., [41,42]). Furthermore, this method was shown to reduce the magnitude of mean errors associated with bathymetric topography estimations by up to 50% (i.e., accuracy of 0.02–0.05 m and precision of 0.06–0.09 m). It was not able, however, to eliminate completely the systematic increase in error with increasing water depth, and difficulties in accurately establishing the position of the water's edge were noted.

In 2017, Dietrich [28] further developed this method by proposing a more sophisticated multi-angle refraction correction approach. Instead of using the orthophotos and DEMs generated by the RPAS–SfM process, Dietrich used the dense three-dimensional (3D) point clouds (which describe the 3D topography of the scene) and applied an iterative refraction correction method based on the positions and view angles of all the cameras used within the photogrammetric reconstruction of the position of each point in the cloud. As found by [27], the challenges of accurately defining the water's edge and the compounding effect on bathymetric accuracy were noted. Dietrich's method requires more inputs than the small angle approach. Guidance on its implementation is available online and is accompanied by a dedicated python script ([43]: pyBathySfM v4.0). Furthermore, Dietrich reports improved error statistics associated with this method of refraction correction, with typical accuracy values of ± 0.01 m and precision values of 0.06–0.08 m. However, a direct comparison of these two approaches using the same input data has not yet been attempted, and nor has either method been tested for quantifying topographic change in the submerged parts of river channels. Therefore, our first objective is to investigate how our RPAS–SfM bathymetric refraction correction methods can be used to accurately measure change in submerged parts of rivers. We include in this a direct comparison of the methods of Woodget et al. [27] and Dietrich [28].

#### *1.2. Water Surface Elevation*

Recognising that one of the key impediments to accurate and reliable refraction correction relates to the challenges of defining accurately the position of the water's edge and the elevation of the water surface, our second objective is to explore different methods for estimating water surface elevations and their influence on the accuracy and precision of below-water topography estimates and estimates of change over time. Knowledge of the water surface elevation is paramount within both refraction correction approaches [27,28], as it forms the lynchpin for computing the apparent water depth, and subsequently the true water depths and refraction-corrected elevations using the geometry of refraction described by Snell's Law.

#### *1.3. Spatially Variable Elevation Error*

The quantification of topographic elevation change in any setting, and particularly within this initial assessment of refraction-corrected bathymetric change, requires the establishment of appropriate error thresholds for each epoch of data and propagation of this error into models of change, in order to distinguish real change from noise [1,16,23,44]. In the past, geomorphic change detection studies have either (a) recognised the limitations of change quantification but not defined or implemented a specific error threshold (e.g., [45]), or (b) have employed spatially uniform error thresholds based on surface topography, data acquisition method, and/or post-processing (known as a 'level of detection' or LoD, e.g., [8,46,47]). However, such broad approaches run the risk of masking real change in some places and detecting change which is not real in others [32,48]. Therefore, as others have done before us, we highlight the need for spatially variable error thresholds [24,32,48–50] and implement them within our work to optimise the accuracy of our change measurements. Previously, spatially variable error thresholds have been based on measures associated with broad polygons of wet or dry areas or of different survey methods (e.g., [8,16,38]), or on spatially-continuous measures of local topographic variability such as roughness, grain size, or slope angle or on sampling point density [1,33,48,51]. However, very few studies have tested the combined impact and significance of multiple proxies on elevation error (with the exception of the fuzzy inference systems used by [1]), and none have yet quantified spatially variable error of RPAS–SfM-derived DEMs in submerged areas. Thus, we take this a step further with our final objective, by utilising the full potential of the exceptionally high resolution RPAS imagery (<0.02 m pixel size) to investigate how this detailed information (specifically concerning topographic complexity, landscape composition, and survey quality/conditions) acquired by each RPAS–SfM survey can help to elucidate the direction and magnitude of elevation errors in submerged and exposed areas. Specifically, we explore how this information might be exploited using artificial intelligence approaches (including multiple regression and a machine learning classifier) to inform the development of an appropriate spatially variable error threshold.

To this end, we study specifically the influence of the following variables on elevation error:

These variables can be obtained for each pixel within the RPAS–SfM-derived datasets and many of them have been shown previously to have some bearing on topographic error [1,26,27,48,49,51]. Elevation error is determined by computing the elevation difference between point cloud elevations and the elevation data acquired from >3500 validation points measured using a combination of GNSS and total station surveys. Given that the existing evidence suggests that simple monotonic relationships between these variables and elevation error is unlikely [1,48], we investigate the use

of multiple regression and machine learning classification algorithms for predicting elevation error. Machine learning (ML) is a form of artificial intelligence (AI), whereby a 'machine' is trained to recognise statistical patterns and trends within data and use them to predict new data outputs. As such, it includes a broad range of algorithms, encompassing everything from simple linear regressions to deep learning-based neural networks [52]. The application of AI/ML for geomorphic error thresholding purposes is novel, with existing studies within freshwater settings applying it predominantly for classification of land cover types from image-derived parameters (e.g., [53–56]) or for the identification of other specific features of interest such as buildings (e.g., [57]) and invasive species (e.g., [58]). As such, our ultimate aim is to create the first high resolution, spatially continuous SfM-derived topographic change models in submerged fluvial environments constrained by spatially variable error estimates.

#### **2. Materials and Methods**

#### *2.1. Site Set-Up*

Our field setting is a 600-m-long reach of the River Teme in Herefordshire, near the England/Wales border, where we collected RPAS imagery in August 2016 and August 2017 (Figure 1). This is a lowland, meandering, dynamic, gravel-bed river system with channel widths of 3–13 m and a total catchment drainage area of 160,000 hectares [59]. In-channel hydromorphic diversity is good, palaeochannels are common (often still connected to the main water course), and median grain size (D50) of exposed surfaces is 21.5 mm. Our section of the River Teme has previously experienced significant channel change following a flood event in the winter of 2014–15 (as evidenced by before and after Google Earth imagery).

Prior to the first RPAS image acquisition survey at the River Teme in August 2016, we established six permanent ground survey markers around the outer boundary of our area of interest (PM1–6, Figure 1). These markers comprised a 50-cm-long steel stake, topped with a flat, plastic yellow survey head and sealing cap which features a plumb point recess for accurate alignment of survey poles. We surveyed the position of these markers in August 2016 using a Trimble R8 RTK GNSS (with a manufacturer-reported accuracy 10 mm in horizontal and 20 mm in the vertical) and collected all subsequent survey data relative to these marker positions. Given the sturdy nature of these markers and light grazing land-use in this area, we were confident that the markers would not have moved between surveys and therefore could be used to ensure the accurate alignment of multiple epochs of data. Despite this, we re-surveyed the markers in August 2017 using a Leica VIVA GS10 GNSS and found that the difference from 2016 marker positions was, on average, less than ±0.008 m, ±0.007 m, and ±0.0185 m in XY and Z, respectively. This difference is very minimal and is within the bounds of uncertainty we would expect as a result of random errors, from pole tilt for example. As a result, we assume no movement has occurred within the intervening time. Three of the six markers could not be found on return to the site in 2017, so an additional two markers were added in new locations and surveyed using the Leica VIVA GS10 GNSS (PM7–8, Figure 1). This allows us to maintain the precision of subsequent internal surveys undertaken using a total station.

**Figure 1.** Location map and site overviews (orthoimagery) from 2017 and 2016 surveys.

#### *2.2. RPAS Surveys*

In August 2016, we flew manually using a rotary-winged DJI Inspire 1 RPAS over the site at a height of ca. 35 m above ground level. High resolution imagery (<0.02 m ground sample distance (GSD)) was acquired using the associated Zenmuse X3 FC350 camera, predominantly at a nadir viewing angle, with a 60–80% level of overlap. The survey was carried out in dry and calm conditions, with sunny intervals. In August 2017, we acquired our second epoch of high-resolution imagery (<0.2 m GSD) using a DJI Phantom 4 Pro RPAS, which was pre-programmed using DJI's Ground Station Pro application, to collect imagery at nadir with an 80% overlap from an altitude of ca. 35 m. Prevailing weather conditions were comparable to those in August 2016. We also collected a number of oblique images from higher altitudes with each survey to incorporate into the SfM processing, to avoid the risk of systematic doming errors which can otherwise result from inadequate camera lens calibration models [27,60–62].

Prior to each survey, we had distributed a number of ground control points (GCPs) throughout the site (*n* = 19 in 2016, *n* = 14 in 2017), ensuring they were positioned to represent the range in elevation distributions across the site (Figure 1). The position of all GCPs was surveyed using a Leica Builder 500 total station in relation to the known positions of the permanent markers.

In 2016, a total of 773 images were acquired, of which 220 images were excluded due to redundancy, blurring, or other visual distortions. In 2017, a total of 1512 images were acquired, of which 471 images were excluded due to redundancy or blurring.

#### *2.3. SfM Processing*

Each epoch of imagery was imported into Agisoft's PhotoScan Pro software for processing using SfM photogrammetry (v.1.4.0 build 5650). Image alignment was undertaken with high accuracy and was optimised following the addition of GCP locations within the model. Exports included a ca. 0.01 m resolution orthophoto, ca. 0.04 m resolution digital elevation model, and a dense 3D point cloud for each survey date. All outputs are referenced to British National Grid (OSGB 1936) using the GCPs and permanent markers established at the site.

#### *2.4. Validation Data*

In both 2016 and 2017, independent topographic validation data were acquired shortly after the RPAS flights in both exposed and submerged parts of the site. Elevation (Z) and XY positions of each validation point were recorded using the total station and tied into the same coordinate system using the permanent markers. Efforts were made to collect validation data using a random distribution, whilst maintaining coverage of all key parts of the area of interest. Particular attention was paid to acquiring validation data over a range of water depths and along breaks of slope, to ensure the validation was not biased to the more easily accessible, flatter parts of the site. This sampling strategy was prompted by knowledge that breaks of slope tend to be areas of greatest error within elevation models [1]. A total of 1522 validation points were collected in 2016, and 2091 points in 2017 (Figure 2).

**Figure 2.** Validation datasets for 2017 and 2016.

#### *2.5. Water Surface Elevations*

Water surface elevations are a critical variable in any refraction correction, since the equations derived from Snell's Law are in terms of water depth rather than absolute elevations. As discussed in [28], any errors in water surface elevations translate directly as an error into the refraction-corrected depths and, ultimately, the refraction corrected elevations. As in other studies, we choose to approximate the water surface as planar in the cross-stream direction (bank-to-bank) and spatially variable in the downstream direction.

For each dataset, we digitised polyline features at the approximate edges of the water surface. The digitised lines tended to be along areas of exposed sediment where the water's edges were clearly visible. In the riffles, where the flows were more complex, the polylines followed the main channel (majority of the flow) and in the side channels, additional point features were digitised at regular intervals (~1–2 m) along the water's edge. The side channel points were not used in the smoothing operation below but were used in the final water surface interpolation.

To extract the elevations along the polyline features, we created a 1 cm buffer polygon on each line segment. The buffer polygons were used to select and extract data points directly from the SfM point clouds in CloudCompare. For the side channels, for all point features, an elevation value was assigned by taking an average of the five nearest point cloud points.

The nature of the extracted point cloud points, with variable elevations at the scale of individual sediment clasts, necessitated that the data be statistically smoothed to generate an average downstream elevation profile that approximated the water's edge. Downstream distances were assigned to all of the extracted points using a stream normal coordinate transformation [63], where the edge points are transformed into coordinates with downstream (along a centreline) and cross-stream (orthogonal to the centreline) components while also retaining their original X, Y, and Z coordinates.

We tested seven different statistical smoothing filters, mostly in the "local polynomial" class of filters. These included Loess and Lowess, along with their robust counterparts, Savitzky–Golay, Gaussian, and Moving Average. The filtering combined both right and left bank points to best estimate the planar cross-stream elevations. All these filters operate via a smoothing window. The optimal window size for each was determined by trial and error, qualitatively judging the fit quality, e.g., preserving known downstream variability, while suppressing high-frequency "noise" and outliers. Once the windows were optimised, the downstream water surface elevation profiles from the seven different filters were compared to determine the best fit for the given set of input data.

To generate the final water surface for each dataset, we resampled the smoothed profiles at 2 m intervals downstream. The 2 m intervals were along the centreline of the stream and we constructed orthogonal cross-section lines that extended beyond the banks at each interval (Figure 3). The attributes of each 2 m interval point were copied to new point features at the both ends of the cross-section lines. We combined the 2 m interval points with the original side channel points and used a Delaunay triangulation to interpolate all the points into a coherent water surface mesh. For the multi-view refraction correction processing, the mesh was used to calculate apparent depths in the point cloud via a Point-to-Mesh distance calculation. For the small angle correction, the mesh was "rasterised" into a 5 cm cell size GeoTIFF, and the apparent depths calculated by raster subtraction. We assessed the performance of this smoothed water surface by direct comparison with a manually digitised water surface (generated as per [27]), through validation of the derived submerged elevations obtained by the small angle refraction correction procedure outlined below.

**Figure 3.** Water surface interpolation points. Centreline points (circles) are generated from the statistical smoothing operations. Edge points (squares) are copies of the centreline shifted orthogonally (black arrows) to create a solid interpolation surface that extends beyond the wetted channel. Additional points are added along the riffles to handle the more complex topography in the side channels.

#### *2.6. Refraction Correction*

In our SfM data, the submerged portions of the point cloud are impacted by the refraction of light at the air–water interface. The result is that the submerged point elevations are too shallow (Figure 4). These points therefore represent the apparent depth (*ha*) based on the water surfaces described above. To correct the point cloud, we used algebraic manipulations of Snell's Law (*n*1 sin *i* = *n*2 sin *r*) to adjust the apparent depth to the "true" depth (*h*).

**Figure 4.** Diagram illustrating through water refraction geometry. Illustrating an off-nadir field of view (FOV) from the camera location (*Xc*,*Yc*,*Zc*) down to the water surface (*WSz*). Point (*Xa*,*Ya*,*Za*) is at the apparent depth (*ha*) and represents the uncorrected SfM point cloud elevations. Point (*Xp*,*Yp*,*Zp*) is the "true" position as the "true" depth (*h*) that the refraction correction is solving. *D* is the point camera distance, *dH* is the elevation difference, *r* is the angle of refraction, *i* is the angle of incidence, *x* is the point/water interface distance, *n*1 and *n*2 are the refractive indicies of air and water. From [28], used with permission.

#### *2.7. Small Angle Refraction Correction*

For the small angle refraction correction process, we imported both the 'old' and the 'new' water surface elevation models into GIS as a raster, alongside the orthophoto and DEM outputs from the SfM process. Following the method of [27], the DEM elevations were subtracted from the water surface elevations to estimate the apparent water depth (*ha*) for each epoch of data. The refraction correction is calculated using the "small angle substitution" by assuming the angles (*r* and *i*) are less than 10◦. From this assumption, the sines in Snell's Law can be substituted with the tangent (sin θ tan θ) d Snell's Law simplifies to the refractive index of water (1.34) multiplied by the apparent depth:

$$h = 1.34 \times h\_a \tag{1}$$

Thus, in GIS, the apparent water depth raster layer is multiplied by 1.34 to create a new, refraction corrected depth raster (*h*). The difference between these two depth rasters (*h* − *ha*) is computed, and the result is subtracted from the original DEM to get corrected absolute elevations. This has the effect of lowering the channel elevations in submerged areas, to counter the underestimation of true bed elevations which occurs as a result of refraction. This process was carried out twice, using first the 'old' water surface and then the 'new' water surface, so that we could compare the two surfaces directly.

#### *2.8. Multi-View Refraction Correction (BathySfM)*

The multi-view refection correction is described in-depth in [28]. With this new research, the software was updated and renamed (now called 'pyBathySfM'), and is available for download on Github [43]. The new version has additional options and a new graphical user interface to make the tool more user-friendly. The minimum inputs for the bathymetric correction are designed to be software agnostic and include: (1) A point cloud dataset (in CSV format, with a water surface elevations attribute added), (2) the camera positions/orientations (X, Y, Z, Pitch, Roll, Yaw), (3) sensor description (focal length, sensor size [physical and pixel size]). The software calculates all of the points that each camera "sees" in its field of view (calculated via a pinhole camera model). For each point–camera combination, the horizontal distance (*D*) and difference in elevation (*dH*) are calculated and used to derive the angle of refraction (*r*). The point elevation (*Za*) is subtracted from the water surface elevation (*WSz*) to get the apparent depth (*ha*). From these variables, the other variables in Figure 4 can be calculated (via basic trigonometric functions) and incorporated into the equations derived from Snell's Law to calculate the true depth (*h*):

$$h = \frac{x}{\tan i} = \frac{h\_d \tan r}{\tan \left[ \sin^{-1} \left( \frac{n\_2}{n\_1} \sin r \right) \right]} \tag{2}$$

where *i* is the angle of incidence, *n*1 is the refractive index of fresh water (1.34), *n*2 is the refractive index of air (1.0). The average true depth for each point can then be subtracted from the water surface elevation to give the corrected bed elevation of the point (*Zp* = *WSz* − *h*).

As discussed in [28], because each point has a fixed apparent depth with variable refraction angles to each camera, there is a different true depth calculated for each refraction angle (camera). The range of values is small at shallow depths and grows with depth. A range of descriptive statistics (e.g., mean, median, mode, percentile ranks) has been explored and we find that the mean value of all of the calculated *h* values for each point provides the best statistic. During the course of this research, we found decreased accuracy when the refraction angles (*r*) were higher (>35◦ off-nadir). These high refraction angles correspond to greater point–camera distances, and also an increase in the amount of water the camera must 'see' through to estimate the bed surface. We decided to incorporate a filtering function into the software to add additional output that uses a limited set of the calculated *h* values, in our case those less than 35◦ off-nadir.

The output of the multi-view refraction correction is a new point cloud (CSV) that retains all the original attributes and new columns for *ha*, *h*, and the corrected absolute elevations (*Zp*).

#### *2.9. Refraction Correction Validation*

Using the validation points that we collected in the field (Figure 2), each of which has precise XYZ coordinates, we assessed the accuracy and precision of our models by subtracting the modelled elevations from the validation point elevations. A small number of outliers were identified and excluded from further analyses. The outliers occurred predominantly in areas of overhanging or dense vegetation over the submerged areas. We then subdivided our datasets into the following categories to quantify the distribution of errors resulting from each approach:

We produced histograms and scatterplots to visualise the error distributions and calculated a variety of error statistics as measures of accuracy and precision. We compared the differences in mean errors and standard deviations between approaches to compare the different water surfaces and different refraction correction approaches. We compared the error statistics for all submerged datasets with those obtained for exposed areas, which to date have typically been much more accurate than those reported for submerged areas [27,28].

#### *2.10. Spatially Variable Errors*

We began our assessment of spatially variable elevation error by extracting several parameters at the location of each validation point, from both the 2016 and 2017 surveys. Table 1 details each of the 11 parameters and how they were extracted from the RPAS–SfM datasets. We chose these variables because they have been used previously for thresholding error when computing geomorphic change (e.g., [1]) or because we suspected that they have a bearing on error based on our prior experience. In addition, we computed the variability of some parameters (slope, roughness, image quality) within a 0.2 m circular window around each validation point, by extracting the maximum, mean, and standard deviation of values within this window. We performed this step using the Focal Statistics tool within ArcGIS (ESRI Inc.). The purpose of this step is to highlight how positioning errors in X and Y translate into elevation (Z) errors, particularly in areas where elevation changes rapidly (i.e., areas of high slope or roughness). For example, in theory, where there is a large range in slope angle and roughness values within the 0.2 m circular window, we expect greater topographic complexity and, therefore, also potentially greater elevation errors.

We tested our theories by assessing individual linear relationships between each variable and the associated error values obtained from our validation data. We also pioneered a new approach for predicting the spatial variability in elevation error, by combining variables together using multiple regression and through the use of machine learning classification algorithms.

We note that, where applicable, we used the DEMs and point clouds derived from the BathySfM–All refraction correction procedure and the smoothed water surfaces to generate the variables specified in Table 1. Given the lack of a consistently superior refraction correction approach (as reported in the results section below), this choice was subjective and one of the other approaches could have been chosen instead.

**Table 1.** Variables used for the prediction of elevation error within the RPAS–SfM surveys. Total number of validation samples for 2016 was 1522 and for 2017 was

#### *2.11. Multiple Regression*

The variables listed in Table 1 were considered for inclusion in a multiple linear regression against the elevation errors computed using the validation data. Variables were assessed for collinearity; where two variables demonstrated a strong correlation (defined as having an R2 greater than 0.7), the variable with the weaker relationship with elevation error (as determined by single linear regression) was excluded from further analyses. This was necessary, given the strong correlations between variables such as roughness and focal roughness. The remaining continuous variables (max. slope focal, min. slope focal, water depth, point density, mean image quality focal, roughness focal, and tie point precision) were standardised by dividing by their standard deviation. Binary variables (vegetation, shadows, reflections) were not standardised. We used MS Excel to run a multiple regression using the chosen variables against (a) elevation errors, where both magnitude and direction of error are considered, and (b) absolute elevation errors, where only magnitude of error is considered. We performed multiple regressions on the 2016 and 2017 validation data separately and combined. We judged the performance of the multiple regressions by considering the multiple R, adjusted R square, standard error, and significance values.

#### *2.12. Machine Learning Classification*

For the first time in geomorphic change detection studies, we also investigated the potential of a machine learning classification approach for predicting elevation error for each validation point using the explanatory variables listed in Table 1. We had initially trialled an ML approach which aimed to predict absolute error values (rather than error categories), but we found that significantly larger training datasets would be required to produce meaningful outputs. Therefore, we binned the values into two sets of error classes: A set of 3 classes (<–0.2, –0.2 to 0.2, and >0.2, representing high magnitude negative error, low magnitude error, and high magnitude positive error) and a set of 10 more detailed classes (<–0.5, –0.5 to –0.2, –0.2 to –0.1, –0.1 to –0.05, –0.05 to 0.0, 0.0 to 0.05, 0.05 to 0.1, 0.1 to 0.2, 0.2 to 0.5, >0.5). This turns the problem into a classification problem, rather than a regression problem.

We started by performing an initial investigation of appropriate classification approaches using the Tree-based Pipeline Optimization Tool (TPOT; [64]), which uses a genetic algorithm to search across multiple common classification algorithms. Manual tweaking based on the TPOT results was used to obtain the best balance of accuracy and pipeline complexity. The resulting classifier was a Gaussian Naïve Bayes classifier using all components of a Principal Components Analysis transform of the explanatory variables as input, using the methods implemented in scikit-learn v.0.20.3 [65]. Separate models were created based on the 2016 validation data, the 2017 validation data, and the combined 2016 and 2017 data. During the ML phase, we observed that the class distributions in the training data were strongly imbalanced (e.g., for the three-class problem with the 2017 data, 90% of the validation points were in the central class). This meant that the classifiers were able to achieve a very high accuracy by simply always predicting outputs within this one class. To combat this, we oversampled the training data using the SMOTE algorithm [66], as implemented in the imblearn v.0.4.3 package [67]. Next, we used a stratified five-fold cross validation procedure for assessing the repeatability of the classification. We observed significant variation in results between runs due to very small numbers in some of the classes (particularly for the ten-class problem). As a result, we ran the five-fold cross-validation procedure five times and averaged the results. After this development and validation of the ML approach, we applied the models to the full raster datasets using the rasterio package [68] to handle raster processing. Finally, we assessed the performance of the ML classifier by comparing predicted outputs with those remaining in the validation dataset to compute overall accuracies and class-specific accuracies, the latter being displayed as confusion matrices.

All machine learning analysis was carried out using Python 3.7.1 using numpy v.1.16.2 [69], pandas v.0.24.1 [70], matplotlib v.3.0.3 [71], and Jupyter Notebook v.5.7.4 [72]. Our code is also available to download [73].

#### *2.13. Error Propagation and Change Detection*

To detect geomorphic change between the two epochs of data, we subtracted the earlier refraction-corrected elevation dataset (2016) from the later one (2017) to create a DEM of difference (DoD). For both epochs, we used elevation data which had been corrected for refraction in submerged areas using the BathySfM–All method. This choice was subjective, given that complex patterns of performance which exist between the different refraction correction methods (Table 5) and another method could have been chosen instead. To isolate real geomorphic change from noise within the derived DoD, we test the implementation of our multiple regression method for predicting spatially variable error. We had initially also intended to threshold the DoDs using the outputs from our machine learning classification approach. However, as detailed in the results, these outputs were not sufficiently robust for the detailed ten-class problem to be used for error propagation purposes. Still, we provide a map of error predictions from the three-class problems for comparison purposes in Figure 9.

First, we used the variables associated with each data epoch (Table 1) to predict the magnitude and direction of elevation error using the multiple linear regressions. We chose to predict both the magnitude and direction of error to maximise the reliability of error propagation during the subsequent quantification of elevation change. We took all input variables as raster layers in GIS and used Equation (3) to predict elevation error (*SVE*) for every pixel in the respective epochs of data, where *MaxSf* and *MinSf* are maximum and minimum slope angles within the 20 cm focal window, *h* is refraction corrected water depth (using BathySfM–All), *D* is point cloud density, *MeanCQf* is mean image quality within the focal window, *MeanRf* is the mean roughness within the focal window, *P* is tie point precision, *V* is vegetation presence, *Sh* is the presence of shadows, *Rc* is the presence of surface water reflections or water surface roughness, σ is the standard deviation of each variable (computed from the validation data for both epochs), *a–g* and *i–k* are the co-efficients associated with each variable, and θ is the intercept.

$$\begin{array}{rcl} SSE &= \mathfrak{a} \Big( \frac{\text{Max} Sf}{\sigma\_{\text{Max} Sf}} \Big) + b \Big( \frac{\text{Min} Sf}{\sigma\_{\text{Min} Sf}} \Big) + c \Big( \frac{h}{\sigma\_{h}} \Big) + d \Big( \frac{D}{\sigma\_{D}} \Big) + c \Big( \frac{\text{MeanCQf}}{\sigma\_{\text{Max}CQf}} \Big) + f \Big( \frac{\text{MeanRf}}{\sigma\_{\text{Max}Rf}} \Big) + g \Big( \frac{p}{\sigma\_{P}} \Big) \\ &+ (i \times V) + (j \times Sh) + (k \times Rc) + \theta \end{array} \tag{3}$$

We performed this computation using the Raster Calculator tool in ArcGIS to output a new raster layer of spatially variable predicted elevation errors for every pixel within each epoch of data (*SVE*2016 and *SVE*2017). We propagated these individual SVEs into the DoD using the 'Error Propagation' tool within the Geomorphic Change Detection v7 plugin for ArcGIS (http://gcd.riverscapes.xyz/). This is based on the approach detailed in Equation (4) [18], where ε*DoD* is the propagated error in the DoD and *SVE*2016 and *SVE*2017 are the spatially variable elevation error predictions for each epoch respectively. Finally, we used the GCD toolbox to threshold the DoD on a cell-by-cell basis to quantify the direction and magnitude of change across the entire site.

$$
\varepsilon\_{\rm DnD} = \sqrt{\left(SVE\_{2016}\right)^2 + \left(SVE\_{2017}\right)^2} \tag{4}
$$

#### **3. Results**

#### *3.1. SfM Modelling*

The georeferencing accuracy and precision statistics associated with each RPAS–SfM survey are presented in Table 2. An evaluation of the validation points indicated that the errors were normally distributed and not affected by any systematic distortions.

**Table 2.** Georeferencing accuracy statistics (all values in metres).

#### *3.2. Water Surfaces*

The statistically smoothed water surfaces, derived from the buffered water's edge points, produced smoothly varying surfaces that provided a better representation of the estimated water surface elevations when assessed visually (Figure 5). This is in contrast to the manually digitised ("point picking") methods used in past studies [27,28]. By picking individual points (from the raster or point cloud) to represent the water surface, there is an increased chance that small variations in the water's edge topography (e.g., topographic roughness of digitising blunders) can have a pronounced effect on the water surface elevations. In Figure 5, we see the effect of these digitising errors as cross-stream lines or also cross-stream mesh facets that translate directly into errors in the calculated apparent depth. These apparent depth errors then get transmitted through the refraction correction to the corrected depths.

**Figure 5.** Comparison of two methods for water surface delineation, manually "picked" edge points (left, e.g., [27,28]) and the method used in this paper using statistically smoothed elevation profiles (right). Colours represent the apparent depth (*ha*) calculated from each water surface method, the black arrows on the manually digitised surface highlight cross-stream errors/artefacts.

For each of the smoothing filters, the window size used for smoothing these data was 151 points (from the buffered point cloud points). Since the downstream densities of points were not uniform, the 151-point window size does not correspond to any particular length scale. For others using this method, we suggest testing several different window scales depending on the stream length, point cloud density, and topographic roughness.

The seven statistical filters varied in their representations of the water surface; the most pronounced differences were at transitions into and out of the steeper sections of the study reach. We calculated the range of water surface elevations that each smoothing filter produced along the interpolated centreline of the river. The maximum difference in 2016 was 0.038 m with a mean difference of 0.007 m. In 2017, the maximum difference was 0.18 m with a mean difference of 0.01 m. To illustrate the effect that the different water surface elevations have away from the centreline, we also calculated the range of values at each of the 2016 validation points. For the calculated apparent depth, the maximum difference between the different water surfaces was 0.074 m with an average of 0.006 m. After refraction correction, the maximum difference in the corrected elevations was 0.033 m with a mean difference of

0.003 m. Overall, a qualitative assessment found that the Savitzky–Golay filter gave the best visual match to the water's edges visible in both years' orthophotographs.

Table 3 presents the quantitative performance indicators for this smoothed water surface (column 3) compared to the manually digitised water surfaces (column 2) used within the small angle refraction correction calculations for both survey epochs. These results demonstrate an unmitigated improvement in accuracy and precision of submerged elevations with the implementation of the smoothed water surfaces. For the 2016 survey, mean error was reduced by ca. 73% and the standard deviation of error was reduced by ca. 69%. For the 2017 survey, mean error was reduced by 50% and the standard deviation was reduced by ca. 25%. Similar improvements can also be observed in the slope angle of the observed versus predicted relationships and the error ranges (as indicated by the min and max errors).

 **3.** Refraction correction accuracy and precision statistics (all values in metres, except slope).

**Table**

#### *3.3. Refraction Corrections*

Table 3 and Figure 6 present the performance statistics associated with the comparison of different refraction correction approaches. The results in columns labelled 3–5 of Table 3 were all computed using the smoothed water surface elevations and Figure 6 presents only the smoothed water surface results. Here, we observe an inconsistent picture of accuracy and precision between the different refraction correction approaches and data epochs. For 2016, the small angle refraction correction produced the highest mean error (0.041 m) and standard deviations of error (0.065 m), indicating a poorer performance overall. However, for 2017, the poorest performance is observed from the BathySfM–All approach. When assessed by the strength of the observed-versus-predicted relationship (i.e., slope of the linear regression line), the small angle correction seems to perform better than the other correction methods in 2016, albeit with the caveat that the regression line has a slightly negative offset (Figure 6). However, this pattern is not evident in the 2017 data. Similar inconsistencies are observed in the minimum and maximum error values. As such, it is difficult to conclude that one particular refraction correction approach is superior for this particular site. It is important to note, however, that almost all of the approaches using the smoothed water surface provide significantly better representations of the submerged topography than the manually digitised water surface. Furthermore, depending on the exact choice of method, these representations are now typically in line with or better than those obtained in exposed areas.

Because the error statistics (Table 3) and distributions (Figure 6) are quite similar, we employed two statistical tests to evaluate if the errors between the different methods were statistically distinct (*HA*) or came from the same parent distribution (*H*∅). The tests were a 2-sample Anderson–Darling (AD) and a 2-sample Kolmogorov–Smirnov (KS), both with a *p*-value threshold of 0.01. For the 2016 corrections, the small angle correction was statistically different than the multi-view corrections using both tests (*p* = 0.001). However, the two multi-view correction error distributions were not statistically distinct using the AD (*p* = 0.017) and KS (*p* = 0.062) tests, meaning that we cannot say with certainty that one method is better than the other. The 2017 corrections were all statistically distinct from one another (*p* = 0.001).

**Figure 6.** Error histograms and scatterplots for 2016 (top half) and 2017 (bottom half) comparing the three methods tested: Small angle correction (left column), unfiltered multi-view correction (BathySfM–All, centre column), filtered multi-view correction (BathySfM–Filt, right column).

#### *3.4. Spatially Variable Error*

#### 3.4.1. Linear Regressions

Table 4 presents the multiple R values for the linear regressions between each variable and each of two measures of elevation error. The strongest correlations are between variables maximum slope within the focal window (0.52), slope (0.50), and standard deviation of slope within the focal window (0.44) with magnitude-only error. Other variables, including water depth, do not show a correlation with magnitude-only elevation error. When magnitude and direction of error are considered, there is a notable lack of correlation between error and any of the variables considered.

**Table 4.** Multiple R values for linear regressions of individual variables with elevation error.

#### 3.4.2. Multiple Regression

Table 5 and Figure 7 present the performance variables associated with the multiple regressions for each epoch of data and epochs combined. All regressions are statistically significant at the 99.9% level. The strongest relationship is between the 2016 variables with magnitude-only elevation error (0.65). Whilst the performance of these multiple regressions (Table 5) is notably better than for the linear regressions for individual variables (Table 4), none of the multiple regression permutations exhibit very strong relationships where the multiple R is greater than 0.7. Standard error is lowest for the 2017 variables and magnitude-only elevation error (0.11 m). Mean residual errors are low (<0.04 m) but the slopes of the observed versus predicted relationships (Figure 7) are far from 1:1. Furthermore, large minimum and maximum residual errors are evident, which push up the standard deviation of residual error (especially for the 2016 dataset). Overall, magnitude-only elevation error is better explained by the variables than is elevation error which has both a magnitude and direction. Combining the 2016 and 2017 datasets degrades the performance metrics in most scenarios.

**Table 5.** Performance of multiple regressions. NC denotes variables which were not computed.

**Figure 7.** Observed versus predicted elevation errors for (**a**) 2016 and (**b**) 2017 datasets, using the multiple linear regression approach for both error magnitude and direction.

Both error magnitude and direction (i.e., positive or negative) are necessary for thresholding quantities of change. Therefore, it was these multiple regression relationships, for the 2016 and 2017 datasets separately, which we used to compute spatially variable error using Equation (3), despite their lower multiple R values. Of the variables included in this predictive equation, we found the presence of vegetation to be the most influential: Its exclusion from the multiple regression reduced the multiple R values from 0.54 and 0.53 to 0.42 and 0.35 for the 2016 and 2017 datasets, respectively. The exclusion of most other variables, one at a time, had minimal effect on the multiple R values for both data epochs, yet the inclusion of only the few most significant variables did not produce higher multiple R values than when all variables were included. As a result, we considered it important to include all variables within the multiple regression at this stage. However, future research would benefit from a more thorough exploration of the most significant combinations of variables.

Table 6 details the standard deviation values used to standardise the continuous numerical variables and the co-efficient and intercept outputs of the multiple regressions (Equation (3)). Figure 8 presents the SVE outputs for each epoch and the resultant propagated error.

**Table 6.** Standard deviations, co-efficients, and intercept values (multiple regression parameters).

**Figure 8.** Predicted spatially variable elevation error maps for 2016 and 2017 datasets (with histograms of residual error associated with these predictions) and propagated error map (subsequently used to threshold elevation change)—all generated using data from the multiple regression approach.

#### *3.5. Machine Learning Classification*

Table 7 presents the accuracy results from the ML classification approach. Relatively high accuracies of ca. 75% and F1 scores are returned for the three-class problem, but significantly poorer values are observed for the ten-class problem. There are small differences in accuracy between the different epochs. More detailed information can be seen in the confusion matrices for the 2017 three-class problem (Figure 9a) and the 2017 ten-class problem (Figure 9b). For the three-class problem, these show good accuracy for the lower two classes (>80% accuracy), but a poorer accuracy for the higher class (with misclassifications split between the other two classes). Figure 10 shows the spatial distribution of error predicted using the three-class scenario. The ten-class problem confusion matrix gives highly mixed results, with reasonable accuracy for the lowest class (–2 to –0.5), middle class (0 to 0.05) and the highest class (0.5 to 6.5). The highest levels of misclassification occur for elevations between –0.1 and 0.1, which are often misclassified as being in the 0 to 0.05 class. The confusion matrices for the 2016 dataset and both years combined show similar patterns.

**Figure 9.** Confusion matrices for the 2017 data classification, for the (**a**) three-class and (**b**) ten-class problem. The matrices are normalised, summing to 1.0 along the rows. Classes are labelled with their lower bound (so in part (**a**), a label of –2 represents the class of –2 to –0.2).

**Figure 10.** Predicted spatially variable elevation error maps for 2016 and 2017 datasets generated using data from the three-class machine learning classification approach.

#### *3.6. Geomorphic Change*

We used the multiple regression approach as a predictor of the direction and magnitude of elevation errors within our two epochs of data, given the low mean errors associated with this approach (Table 5). We chose not to use the outputs of the ML classification approach on this occasion, given the poor accuracies associated with the ten-class outputs. Whilst the classification results for three-classes are much better, our aim was to predict error with greater precision than is permitted by only three categories of error.

Figure 11 presents a map of geomorphic change at this section of the River Teme, thresholded by our error estimates using the multiple regression approach. Associated statistics which quantify areal and volume change are presented in Figure 12. Over the 1 year period between our surveys, this section of the River Teme experienced a net volume increase of 1155 m3, as indicated by the teal-coloured areas within Figure 11. We note that a proportion of this increase will relate to growth of vegetation (as seen in the vicinity of trees when comparing Figure 11 with Figure 1) and does not solely represent sediment deposition.

**Figure 11.** Elevation change map (2016–17), thresholded by elevation errors predicted using the multiple regression approach.

**Figure 12.** Quantities of elevation change (2016–17) in terms of volume (**left**) and area (**right**). Note that the x-axes have been truncated for clarity purposes. Total areal lowering was 8097 m2, total areal raising was 17,570 m2, and net volume difference was +1155 m3.

#### **4. Discussion**

#### *4.1. Refraction Correction*

For the first time, our research has demonstrated that our RPAS–SfM refraction correction methods are capable of quantifying above- and below-water fluvial geomorphic change at hyperspatial resolutions, without the need for different methods between exposed and submerged areas or for acquisition of calibration data in submerged areas. This represents a significant time and cost saving compared to existing approaches (e.g., [37]) and has positive implications for surveyor safety, as there is no need to enter the water.

Our accuracy statistics (mean error) indicate that the below-water topographic reconstructions are now broadly as accurate as those obtained in exposed areas, and our precision values (standard deviation) demonstrate that they are significantly more reliable than for exposed areas (Table 3). This may result from a lack of vegetation in submerged areas which generates noise within the point cloud, thereby reducing precision of our topographic reconstructions in much of the exposed parts of the site.

Our results show no consistent pattern of compromise on the accuracy or precision of the topographic data compared with approaches which use different methods in submerged areas and/or those which require calibration data. For example, the accuracy and precision values we report here compare favourably with those reported by [37,42] (Table 8), where optical calibration approaches or site-specific refraction correction coefficients are used in submerged areas—both of which require extensive calibration data collected in-channel. Comparison of our results with those of [38] present a mixed picture, with our mean error and standard deviation values being higher but RMSE values being broadly comparable. Further exploration of this is restricted by the limited density of validation points and a lack of information on validation point distribution and water depths within the work of [38].

**Table 8.** Comparison of RPAS-based approaches to below-water topographic reconstruction and change detection.

For the first time, our research has also demonstrated minimal differences in performance of the small angle and multi-view refraction correction methods applied at this site [27,28]. This suggests that either of these methods may be implemented by those wishing to produce similar outputs at their sites, assuming that the basic caveats for successful implementation are met (i.e., predominate camera angles [nadir vs. off-nadir], clear water, water depths < ca. 1.5 m, minimal water surface roughness). Those with a background in GIS may prefer to implement the small angle correction, whilst those familiar with running Python code may prefer the multi-view approach (the latest version has a user-friendly graphical interface). Here, we prefer to use the multi-angle approach as detailed in our GitHub repository [43] as it offers the flexibility to work with large, high-resolution 3D point clouds (rather than 2D rasters) and modify by parameters such as camera view angles, as we used within our BathySfM–Filtered scenario. We suggest that future research should focus on the comparison of these methods over a wider range of study conditions. In particular, establishing, quantitatively, the range of settings in which our refraction correction methods can be applied successfully is of importance. Dedicated laboratory and field experiments which aim to determine the limits of water depth, clarity, and surface roughness would be beneficial and would allow a quantification of the percentage of river environments in which our methods are applicable. Such data will be key for understanding the value of our methods for routine river monitoring and management purposes on a larger scale. Further exploration of the importance of camera view angles would also help inform routine RPAS surveys in such scenarios.

#### *4.2. Water Surface Elevations*

Our second key finding is that accurate water surface elevations are critical for accurate refraction correction in submerged areas and are capable of reducing mean elevation errors by 50–73%. Therefore, we recommend the use of our improved water surface estimation approach for all those wishing to reconstruct submerged fluvial topography using RPAS–SfM derived data.

Specifically, our investigations suggest the use of a Savitzky–Golay filter within a smoothing window produces the most accurate results; however, the low variability between the smoothed water surfaces would suggest that the local polynomial regression filters in general (Savitzky–Golay and normal/robust Loess and Lowess filters) are a good choice for characterising the water surfaces. The other filters that we tested (Moving Average and Gaussian) accounted for a large portion of the variability reported in the results, so these filters are not recommended.

One of the key variables in developing these filters is the smoothing/sampling window. In our dataset, the best results came from a 151-point smoothing window. We settled on this window via a trial and error process, along with a qualitative assessment of the "best-fit" of multiple smoothing window sizes. For other datasets, we would encourage researchers to test multiple filter windows to find the most appropriate size for their data. One of the challenges in choosing the correct window size was that smaller window sizes had a tendency to create "uphill" water surfaces in certain sections, influenced by small outliers in the edge point data. In our case, the larger window sizes helped moderate the influence of small-scale variations in the edge points.

Whilst this methodology represents a notable improvement over past studies (i.e., [27,28]), there remain some limitations. The method still relies on having clearly identifiable banks and bars to have sufficient data to perform the statistical smoothing. In river environments where the banks are heavily vegetated or the channel is incised and there are no visible banks, extracting enough information to utilise this method will be challenging. The assumption of a cross-stream planar water surface is a reasonable approximation for most streams, but there are situations where deviations from this assumption could impact the results. Super elevation of the water surface on the outside of meander bends would cause water surface elevations to be higher on the outside of the bend and potentially lower on the inside. In-stream features that cause local hydraulic disturbances (e.g., exposed and submerged sediment clasts, standing waves, or dune bedforms) will also cause locally variable water surface elevations that are not accounted for in our approximation.

Future work should focus on validating this methodology with detailed surveys of the water's edges, along with direct measurement of the water surface. The ultimate goal will be to develop methods and/or sensors that have application across a wide range of river types and that can be used to extract the water surface elevations, thereby enabling researchers to apply bathymetric corrections over longer river segments without the need for in situ water surface measurements. This will subsequently permit the collection of extensive measurements of river system variables.

#### *4.3. Spatially Variable Error for Geomorphic Change Detection*

Our third key finding is a lack of linear correlation between any individual variable and elevation error. This finding is not surprising, given numerous factors are known to contribute to such errors. It reinforces the existing observations that levels of detection used to threshold geomorphic change (in both magnitude and direction) in such settings should not be informed by single variables alone [1,48]. Furthermore, the lack of correlation between error and water depth evidences the success of our improved water surfaces, combined with our refraction correction approach. Previous studies had failed to break down this correlation so successfully (e.g., [27]).

Fourth, we find that multiple linear regressions (which regress elevation error against a combined total of 11 RPAS–SfM-derived variables) are capable of explaining up to 54% of the magnitude and direction of elevation error in our topographic surveys. This is in a scenario where we have controlled for sources of error within the SfM process in so far as is possible/known at the current time (including use of sufficient GCPs which are adequately distributed, a precise internal control network, and both oblique and nadir RPAS imagery to mitigate against systematic camera lens distortions). These multiple regressions allow us to model the spatial variability of error with mean accuracies of 0.01–0.04 m and standard deviations of 0.14–0.22 m. To our knowledge, this is the first such model of spatially variable elevation errors in RPAS–SfM surveys in a geomorphological setting, which allows us to quantify and threshold geomorphic change with a known level of accuracy and precision.

In the last few years, the RPAS–SfM approach has been described as a 'game-changer' for quantitative 3D analyses in geomorphology. Despite this, comprehensive, spatially-continuous, quantitative estimates of associated elevations errors at very high spatial resolutions have been lacking. Thresholding calculations of geomorphic change using spatially variable elevation errors is critical for the accuracy and reliability of subsequent quantitative analyses and decision making in a management context. Therefore, we anticipate that with further testing and development of our model over a range of different topographic settings, landscape compositions, and survey qualities/conditions, it has the potential for routine application for error estimation within the RPAS–SfM workflow. However, further testing and development is required to improve on the current limitations of our approach. For instance, we note that 46% of elevation error remains unexplained by our model, and therefore unaccounted for within our maps and statistics of geomorphic change. Whilst mean residual errors are low, we also report notably large minimum and maximum residual elevation errors within our resulting predictions of spatially variable error (Table 5). Such large residual errors are rare (e.g., 85% (2016) to 92% (2017) of validation points have residual errors of ± 0.2 m) but deserve greater attention in future, as they tend to occur in areas of greater geomorphic complexity. Furthermore, we observe a significant lack of a 1:1 relationship between observed and predicted elevation error, as indicated by the slope values in Table 5 and Figure 7. We suspect that these poorer results relate to a number of key limitations to our approach:

Our method assumes a linear relationship between error and the included variables. Whilst this may be appropriate for some variables, such as slope or roughness, it makes less sense for the binary variables, such as vegetation and reflection/surface roughness. In future, rather than noting the presence or absence of vegetation, it may be worth exploring the inclusion of a continuous numerical variable like vegetation height above ground level. Whilst this is less easily computed using RPAS–SfM surveys than airborne LiDAR, for instance, the interpolation or extrapolation of terrain elevations across vegetated areas may go some of the way to ameliorating this current issue. This is likely to be easier in areas of reduced geomorphic complexity, such as for relatively flat floodplains rather than deep-seated landslides.

Our method does not account specifically for elevation error resulting from uncertainties in the co-alignment of the RPAS–SfM surveys with the calibration/validation data. This gives rise to the well-documented situation whereby the magnitude of elevation error is strongly influenced by slope (where relatively small offsets in X and Y result in large errors in Z). However, slope does not lead to a consistent positive or negative elevation error because the orientation of the slope varies in relation to the direction of the XY offset. This goes part of the way to explaining why our multiple regression, which predicted magnitude-only error, outperformed that which attempted to predict both magnitude and direction of error. Milan et al. [49] and Carbonneau et al. [74] suggest that the alignment of elevation data as point clouds using cloud-to-cloud registration or the M3C2 algorithm [21] within CloudCompare can ameliorate this issue and suggest that in future, the calibration/validation data used to develop and test models of spatially variable error should be available as point clouds. We anticipate that calibration/validation data which is larger and denser than our total station survey may be required here, such as that obtained using terrestrial laser scanning. This should be a priority for subsequent development of our model.

To give an indication of the impact of these limitations, we re-ran our model, this time excluding all validation points which fell within areas of dense vegetation or water surface reflections/roughness and/or on slopes of greater than or equal to 70 degrees. The output was a reduction in mean residual error of elevation error predictions from –0.04 to –0.03 m for the 2016 dataset and from –0.01 to –0.003 m for the 2017 dataset. Improved precision values were also observed: from 0.22 to 0.11 m (2016) and 0.16 to 0.09 m (2017). The slopes of the observed versus predicted elevation error relationships (Figure 7) were reduced from 0.1297 to 0.0327 (2016) and from 0.1666 to 0.0964 (2017). This suggests that, as we begin to remove the larger elevation errors from the model (resulting from dense vegetation, steep slopes (i.e., XY offsets) and reflections/roughness), we remove the systematic offsets which the model has 'learnt' and are left with what is likely to be random noise.

Whilst we can realistically aim to reconstruct systematic errors (or bias) within models, we will never be able to predict random errors (noise). For this reason, [1] argued that a deterministic model of elevation error can never be constructed unambiguously. As a result, they use a fuzzy inference system for error estimation based on membership of their GPS validation data to three functions: Slope, point density, and point quality. We suggest that a deterministic prediction of error is possible to an extent and, so long as we keep in mind the occurrence of random error in addition to the error we predict, that implementation of a model such as ours is appropriate. With the advent of artificial

intelligence, we believe it is only a matter of time and effort, in collating sufficient and appropriate evidence, before machine/deep learning approaches are capable of predicting systematic errors within RPAS–SfM elevation datasets routinely.

Our initial findings suggest that variables associated with topographic complexity and landscape composition (e.g., slope, vegetation) have a greater bearing on elevation error than do variables associated with the SfM process itself (e.g., tie point precision). However, we also note that the relationships between these different variables are complex. For example, the SfM process relies on image texture for tie point identification, yet texture is itself determined by the composition and complexity of the landscape being surveyed. James et al. [49] argues that point cloud topographic datasets, thresholded by SfM tie point precision and aligned using the M3C2 approach produce "the most plausible distribution of topographic change" ([49], p. 1782). These results could be interpreted to mean that SfM tie point precision is the most reliable predictor of elevation error in RPAS–SfM surveys. However, whilst [49] accounts for horizontal offsets in this approach, they do not specifically explore the influence of other variables on elevation error. Furthermore, their judgement of the most 'plausible' approach is based on qualitative expectations of change at their site, rather than quantitative evidence collected contemporaneously. As such, we suggest that further investigation into the relative importance of a wide range of different variables on elevation error is required, with particular focus on how we might reduce or eliminate such errors. Within this paper, we have demonstrated such an approach for reducing topographic error which results from the presence and depth of water.

Our final key finding is that the machine learning classification approach implemented here is successful only in predicting elevation error within three broad classes: 'High negative error', 'high positive error', and 'low error', with overall accuracies of >75%. Whilst such a classification may be helpful in some broader scale analyses of change, our aim was to develop a method of quantifying spatially variable elevation error more precisely in order to threshold measures of geomorphic change at the high spatial resolutions provided by RPAS–SfM surveys. We anticipate our approach was significantly restricted by the size of the training dataset, which was particularly lacking in data points representing the larger elevation errors. Within the geomorphic-remote sensing community, the collection of validation/training data using traditional surveying approaches (e.g., total station/GPS surveys) is well established. However, as we move towards the exploration of more complex artificial intelligence (AI) tools for developing such predictive relationships, it is likely that we will need to acquire validation/training data on a completely different scale. That is, AI/ML approaches typically require millions of data points to 'learn' from, rather than the hundreds or thousands we might employ for standard regressions. Acquisition of such large datasets at scales comparable to high resolution RPAS surveys will require the use of terrestrial or airborne laser scanning. We anticipate that a more advanced ML approach, of which we are already seeing more widespread application within geomorphological fields (e.g., [75]), in combination with a significantly larger training dataset, will facilitate the development of a much improved AI model for predicting spatially variable elevation error for any given topographic dataset derived from the RPAS–SfM workflow.

#### **5. Conclusions**

In this research, we have shown that we are able to quantify submerged geomorphic change at levels of accuracy commensurate with the exposed areas of our river system with a single photogrammetric approach and without the need for calibration data. For the aerial imagery we collected (predominantly nadir views), there are minimal differences in results produced by different refraction correction procedures (small angle vs. multi-view), which may allow users a choice of software packages/processing complexity depending on their research objectives. Our improvements to estimations of water surface elevations with statistical smoothing in the downstream direction are critical for accurate topographic estimation in submerged areas and can reduce mean refraction corrected elevation errors by up to 73%. We have shown that machine learning, in the form of multiple linear regressions between error and 11 independent variables, can be used to generate high resolution, spatially continuous models of geomorphic change in submerged areas, constrained by spatially variable error estimates. Our machine learning model is capable of explaining up to 54% of magnitude and direction of topographic error, with accuracies of less than 0.04 m. With on-going testing and improvements, this machine learning approach has potential for routine application in spatially variable error estimation within the RPAS–SfM workflow.

**Author Contributions:** Conceptualization, A.S.W.; formal analysis, A.S.W., J.T.D. and R.T.W.; funding acquisition, A.S.W.; investigation, A.S.W. and J.T.D.; methodology, A.S.W. and J.T.D.; software, J.T.D. and R.T.W.; supervision, A.S.W.; validation, A.S.W. and J.T.D.; writing—original draft, A.S.W., J.T.D. and R.T.W.; writing—review and editing, A.S.W. and J.T.D.

**Funding:** This project was supported in part by funding from the University of Worcester (UK), where Amy Woodget was previously based.

**Acknowledgments:** We gratefully acknowledge the field assistance of Paul Attwood, Robbie Austrums, James Atkins, Fleur Visser, Mark, Sally and William Woodget, Carlos Puig-Mengual, Cat Fyffe, Steve Dugdale, Katie Long, and John Powell. We thank Mike James, of Lancaster University, for computing tie-point precision of our SfM surveys using his 2017 approach, and Olivia Wilson, of the Ordnance Survey, for expert guidance with the machine learning classification problem. Constructive comments from Patrice Carbonneau, Leonardo Camelo, and three anonymous reviewers helped to improve an earlier version of this manuscript. We also thank the landowners. The code produced in this publication is available online as follows. Machine learning: https://github.com/robintw/ WoodgetDietrichWilson2019, Multi-view refraction correction: https://github.com/geojames/pyBathySfM. Please note that we do not provide the RPAS imagery associated with this article, due to landowner privacy concerns.

**Conflicts of Interest:** The authors declare no conflicts 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/).

## *Review* **Remote Sensing of River Discharge: A Review and a Framing for the Discipline**

**Colin J Gleason 1,\* and Michael T Durand 2**

Received: 13 February 2020; Accepted: 25 March 2020; Published: 31 March 2020