Potential of a Non‐linear Full‐Waveform Stacking Technique in Airborne LiDAR Bathymetry

Airborne LiDAR bathymetry is an efficient measurement method for area-wide acquisition of water bottom topography in shallow water areas. However, the method has a limited penetration depth into water bodies due to water turbidity. This affects the accuracy and reliability of the determination of water bottom points in waters with high turbidity or larger water depths. Furthermore, the coverage of the water bottom topography is also limited. In this contribution, advanced processing methods are presented with the goal of increasing the evaluable water depth, resulting in an improved coverage of the water bottom by measurement points. The methodology moves away from isolated evaluation of individual signals to a determination of water bottom echoes, taking into account information from closely adjacent measurements, assuming that these have similar or correlated characteristics. The basic idea of the new processing approach is the combination of closely adjacent full-waveform data using full-waveform stacking techniques. In contrast to established waveform stacking techniques, we do not apply averaging, which entails low-pass filtering effects, but a modified majority voting technique. This has the effect of amplification of repeating weak characteristics and an improvement of the signal-noise-ratio. As a consequence, it is possible to detect water bottom points that cannot be detected by standard methods. The results confirm an increased penetration water depth by about 27% with a high reliability of the additionally extracted water bottom points along with a larger coverage of the water bottom topography.


Introduction
Water bodies are often subject to significant dynamic changes and may therefore require continuous monitoring. On the one hand, changes in water bottom topography are the result of human activities, e.g., river regulation or cross structures for shipping to raise water levels (Anderer et al. 2011). On the other hand, riverbeds are subject to natural changes due to erosion and sedimentation caused by the flow of water, river floods, or sudden changes caused by natural disasters (Charlton 2007;Benito and Hudson 2010). Thus, permanent monitoring of the riverbed is a required waterway management task to guarantee the good condition of the river. If necessary, hydraulic engineering interventions such as lowering of groins, clearing out groin fields, or widening of streaklines might become necessary. Similar effects occur in coastal areas after storm or flood events. For coastal protection purposes, monitoring of these morphologic changes is important (Christiansen 2016). Continuous measurement or monitoring of the water bottom topography is especially necessary in shallow water areas of the heavily trafficked North and Baltic Seas to guarantee safe navigation for shipping.
In Europe, the European Water Framework Directive (EuropeanUnion 2000) specifies that inland and coastal waters (surface waters) are to be protected, enhanced and restored [Art. 4, "Environmental objectives" in Europe-anUnion (2000)]. In addition, the EU member states are responsible for the permanent monitoring of the condition of the waters. For inland waters, for example, this includes the aforementioned active waterway management. Regarding the maritime areas, topographic surveying of the seabed and tidal flats, recording of underwater obstacles, and mapping of the coastline and other objects in the coastal area play an important role.
The most common and established measuring method for water bottom acquisition is shipborne multi-beam echo sounding (MBES) (Lurton 2002). Georeferencing is done by high-precision GNSS and IMU systems. The acquisition of large areas is very time consuming and therefore expensive, but due to its penetration depth in deeper water sections, it is indispensable for hydrographic surveys. Airborne LiDAR bathymetry (ALB) represents an alternative for the acquisition of water bottom topography of shallow water areas, where echo sounding is rather inefficient. The water body is scanned by ALB systems with green laser pulses. To support the water surface detection, optional use of a near-infrared laser is possible (Guenther et al. 2000). The green laser pulse is emitted from the ALB system at a known position and interacts with all objects along its path such as the water surface, small particles in the water column, and the water bottom ( Fig. 1). In addition, the water itself affects the propagation speed of the laser beam. The backscattered signals arriving at the receiver are recorded by the ALB system's detector. The position and orientation of the laser scanner as well as the deflection of the emitted laser beam are also provided by high-precision GNSS, IMU sensors, and the mirror angle.
ALB systems come with the advantages of extensive and time-efficient water body measurement. However, the method has a limited penetration depth into the water body, affected by interactions of the emitted laser pulse with the water column. The signal strength is weakened by water turbidity to a high degree (Guenther and Goodman 1978). Diffuse scattering and absorption attenuate the signal strength within the water column (Guenther 1985). Furthermore, the scattering characteristics and the reflectivity of the water bottom lead to an additional attenuation of the signal. Finally, the power level and the divergence of the laser beam are device-specific properties that also have an influence on the signal energy. All of these facts have a great impact on full-waveform data processing and the subsequent detection of water bottom points. Figure 2 illustrates the issue using two exemplary signals. Figure 2a represents the ideal case with distinctive signal echoes from the water surface and water bottom enabling reliable detection of the water bottom echo. In contrast, Fig. 2b shows a typical full-waveform in turbid or deeper water bodies without a pronounced ground echo, where the bottom is thus not always reliably detectable.
In Richter et al. (2017), we presented an approach to improve the bottom echo detection in deep water areas by correcting the laser pulse history for attenuation effects. The approach is based on the individual processing of the recorded full-waveform signals. First, the signal decay coefficient is determined from the individual full-waveforms.
Based on these coefficients, a correction model is applied on the raw full-waveforms, removing the turbidity effects from the received signal. As a result, weak water bottom echoes are amplified, facilitating the detection of water bottom points. Obviously, the amplification affects both ground echoes and noise, and thus, the developed method is limited to data sets with a high signal-to-noise ratio. In the presence of stronger noise, the amplification of the noise components in the signal again prevents reliable detection of weak bottom echoes.
Given that in an isolated analysis, a weak ground echo is hardly or not at all distinguishable from the detector's noise level, a reliable detection of these weak water bottom echoes requires further information. The combination of full-waveform data from closely adjacent measurements in the form of full-waveform stacking is a valuable possibility. In the literature, there are already various applications for the stacking of signals for peak detection purposes. Stilla et al. (2007) used full-waveform stacking in terrestrial laser scanning for weak pulse detection, which is below a noise-dependent threshold of pulse detection methods. In Plenkers et al. (2013), the stacking of seismic signals serves to detect local weak seismic activities. Roncat and Mandlburger (2016) use full-waveform stacking for the detection of water surfaces because the detection does not work reliably for data with a low signal-to-noise-ratio when using single laser shots. All of these applications clearly show the advantages of combining measured signals, namely, (1) the detection of weak events in recorded signals and (2) the increase or improvement of the signal-to-noise ratio. In Mader et al. (2019), the successful application of full-waveform stacking techniques to ALB measurement for water bottom detection and extraction was presented. By stacking closely adjacent full-waveforms, it was shown that the signal of weak water bottom peaks could be amplified and at the same time, the noise signal could be reduced. The methodology was applied iteratively from the shore to the middle of the water body. However, this approach was not optimal due to the unreliable determination of very shallow water depths that could lead to errors in further water bottom determinations.
This article presents an advancement of the approach presented in Mader et al. (2019) for processing bathymetric fullwaveform data considering the correlated neighbour information and its application to a part of the Elbe River with a thorough evaluation of the results. Various methods were embedded in an extended processing scheme with the goal of improving the results regarding the detection of deeper water bottom points and the consequent increase of water bottom coverage. For this purpose, closely adjacent measurement data (full-waveforms) are evaluated and combined into a stacked full-waveform. Weak ground echo signals are amplified due their redundancy, and at the same time, the noise influence on the signal is reduced, allowing a more reliable detection of weak ground signals. This information is used in the individual full-waveform measurements for the extraction of the bottom signal and the associated determination of the water bottom point. Furthermore, analysis of the achieved results regarding the deviations compared to reference data (e.g., echo sounder measurements) is also an important aspect to prove the overall reliability of the approach.
The remainder of the paper is structured as follows. Section 2 contains information on the survey area, the sensor specifications, and the data properties. In addition, the reference data and their accuracy for the evaluation of the results are presented. Section 3 explicates the applied methods with the main focus on combining the measurement data to realize the use of correlated neighbourhood information, i.e., the full-waveform stacking. Furthermore, this section details the tools for analysis, evaluation and filtering of the full-waveform stacking results, which should guarantee a high reliability. Section 4 presents the results of the extended processing algorithm using mean deviations, root mean square, mean absolute deviation (Hampel 1974;Sachs 1982), median absolute deviation (Sachs 1982) and reliability information compared to reference data. The article ends with a discussion about improvements and limits of the approach and a conclusion of this contribution.

Study Area and Data
The development and assessment of advanced full-waveform processing methods requires appropriate high-resolution full-waveform measurement data from a water body with adequate turbidity regarding the expected water depth. This means that a high expected water depth requires a low turbidity, and a low water depth allows a higher turbidity. Ideally, the evaluation of the developed methods is based on reference data for the water bottom topography that have superior accuracy. In our study, we used a data set of the Elbe River provided by the German Federal Waterways and Shipping Office Dresden (WSA) and the German Federal Institute of Hydrology (BfG). The study area is located between the German towns of Klöden and Elster. The data did not completely fulfil the criteria mentioned above, as highly accurate reference data are usually not available in natural water bodies. Nevertheless, the provided data have a high spatial resolution, and the water turbidity was sufficiently low at the time of data acquisition.
In addition to echo sounder measurements, the results of the standard processing method (online waveform processing) serve as validation data for evaluation purposes. The points of the online waveform processing contain the water bottom that reaches up to a water depth of approximately 1.65 m. Overall, the available data were well suited for the development and evaluation of our new processing methods. In the following, the study area, the ALB measurement system, the characteristics of the measurement data, and the comparative data used for control purposes are presented.

Study Area and Measurement Data
The ALB data acquisition of the study area took place in spring 2015 using a RIEGL VQ-880-G (RIEGL 2018) with a point density ranging from 11 points/m 2 up to 55 points/m 2 in a flying height of approximately 380 m above ground using a Palmer scan pattern with a 20 • incidence angle. In addition to the 3D point cloud, the full-waveform was digitized and recorded (RIEGL 2017). Each full-waveform consists of 60-200 samples with a constant sample time interval of 0.575 ns that corresponds to an object distance under water of approximately 6.5 cm for a light velocity in water of 225,000 km s −1 .
For the development of the extended processing approach, several subregions were selected (Fig. 3). In this contribution, evaluation results of Elbe 3, 4, and 5 are presented. The properties and results of the areas Elbe 1, 2, and 3 are very similar, so the results of Elbe 3 are representative for these areas. In addition, Elbe 3 features a large number of validation points from ALB online processing and thus offers good control conditions. The Elbe 4 area includes deeper Digitized ALB waveforms with a two clearly detectable echoes for water surface and water bottom. In b, water surface echo is clearly visible, but the weak water bottom echo is not reliably detectable subareas, and for Elbe 5 as well as Elbe 6, the measurement density is clearly lower (see Table 1) than the other investigation areas so that the results of these subareas allow conclusions about the potential and reliability of the new processing approach under rather difficult conditions.

Reference Data
To assess the results of the processed data, two data sets are available: the results of ALB online processing and echo sounder measurements. The echo sounder data were acquired in autumn 2015 by a vessel. This measurement vessel is equipped with a positioning system (Leica GX1200 RTK SAPOS Rover), an orientation system (Seatex MRU5), and an electroacoustic sounding frame (status: 2015). In detail, the electroacoustic frame consists of 37 single-beam transducer, which are distributed over a span of 12 m on two cantilevers. The echo sounders (manufactured by Dr. Fahrentholz, Kiel) operate with a frequency of 200 kHz and exhibit a depth resolution of 1 cm similar to the Kongsberg EA MCU system (Kongsberg 2009). The 3D point accuracy is stated in Weiß and Wirth (2015) with approximately 7.5 cm in embankment areas and 6.4 cm in the areas of the flat riverbed.

Methods
The basic idea of our approach is the use of spatially adjacent information for a reliable detection of weak water bottom echoes in full-waveform data. The measured data are processed in several processing steps (Fig. 4), starting with the partitioning into processing units. This defines the immediate neighbourhood, which depends on the data density and the ground topography of the water body (Sect. 3.1). Beyond that, prior water-land classification of the processing unit reduces the processing effort on the water body units.
Following the partitioning, the full-waveform data are processed and the extraction of water bottom points from each individual full-waveform is performed. The core method of full-waveform processing is the full-waveform stacking approach, i.e., the combination of all full-waveforms of a processing unit to one stacked full-waveform. Each stacked full-waveform is analysed, and a search corridor for water bottom echoes is estimated. As a result, each processing unit has exactly one search corridor for further processing. These search corridors are applied to the individual full-waveforms to detect the water bottom echoes and calculate the water bottom points. A search corridor is a small section in a full-waveform in which water bottom peaks are searched. It is defined by the number of samples between water surface peak and water bottom peak (number of samples) as well as by the width of the water bottom peak (peak width, Sect. 3.4) as a tolerance range (Fig. 5). The full-waveform processing is performed for every processing unit (Sect. 3.1) and is divided into the following steps: (1) detection of water surface echoes in all individual fullwaveforms using analysis tools: detection of all local maxima (Sect. 3.2) and determination of their significance (Sect. 3.3);

Fig. 3
The study area (part of the Elbe River) including subregions for development and evaluation of the extended approach. The results of the evaluation are shown for the red-coloured areas in this contribution. The grey-coloured areas were evaluated, but the results are similar to Elbe 3, 4 or 5 and are therefore not shown in this contribution. Blue areas present the acquired water bodies (Map data ©2020 Google) (2) full-waveform stacking over all full-waveforms of a processing unit (Sect. 3.5); (3) analysis of the stacked full-waveform: detection of all local maxima and determination of their significance; (4) determination of the bottom echo in the stacked fullwaveform and derivation of the number of samples as well as peak width (Sect. 3.4); (5) filtering of full-waveform stacking results including directly adjacent surrounding processing units using the search corridor as a check corridor (Sect. 3.6); and (6) detection and extraction of the water bottom echo as well as the corresponding water bottom point in every individual full-waveform within the estimated search corridors of the stacked full-waveforms (Sect. 3.7).
All individual work steps are described in detail below. In the following, the water depth of a full-waveform is equal to the number of samples between the water surface and water bottom echo considering the laser beam passing obliquely the water column, the propagation speed of the light under water, and the sample time (= time from sample i to sample i+1).

Partitioning of Survey Data
The detection of weak water bottom echoes assumes that the measurement data are closely adjacent and have similar characteristics. For this, a suitable definition of neighbourhood is necessary. Based on the available data, the following partitioning methods could be considered, which refer only to the points of the water surface (regardless of the corresponding water bottom point's position): (1) use of all data within a maximum distance in the lateral direction; (2) use of surrounding vertices of a certain feature point of a triangulated irregular network (TIN); (3) partitioning in a regular grid; and (4) a flexible adaptation of the grid cell size depending on water depth and turbidity.
The abovementioned possibilities differ in their applicability with regard to point distribution, point density, characteristics of the water bottom, and calculation effort. The first two cases are applicable to all kinds of point distributions but entail a high processing effort. For each full-waveform, the neighbouring full-waveforms must be searched, and the full-waveform stacking is carried out. The high processing effort is caused by calculation of a stacked full-waveform for each individual full-waveform. The application of a TIN ensures the use of closely adjacent information, which is well distributed about the whole horizon. In the third case, a uniform point distribution and density are required. The full-waveforms must be assigned to the regular grid, and subsequently, one stacked full-waveform per grid cell will be estimated. As a consequence, the processing effort is low.
In the fourth case, the grid cell size is adapted locally to the uneven point distribution. The processing effort increases because the point distribution and density must be analysed for the entire area beforehand. The present data set of the Elbe River has a rather uniform point distribution and density, so utilization of a regular grid is sufficient. The non-uniformity of the point distribution of a Palmer scan pattern is taken into account by choosing a sufficiently large grid cell width. An important parameter for using a regular grid is the grid cell size, which depends on the data density and the ground topography characteristics (strong or quick changes in terrain slope). For optimal processing results, the following facts should be taken into account when choosing the grid cell size: -a too small grid cell leads to a low redundancy of fullwaveform data, which increases the risk of incorrect water bottom detections; -an extension of the grid cell size results in a higher number of full-waveforms and therefore a higher reliability of the results; and -a too large grid cell can affect a weakening of the water bottom echo in the stacked full-waveform because the correlation or similarity of full-waveform data decreases with increasing spatial distance. As consequence, water bottom echoes may not be detected.
An optimal grid cell size was empirically determined by testing different grid cell sizes and comparing the number, accuracy, and reliability of the detected points. The results are shown in Sect. 4.1.

Full-Waveform Analysis-Peak Detection
There are different approaches for the detection of backscattered object signals in airborne laser scanning. In Wagner et al. (2004) and Wang et al. (2015), several conventional peak detection methods are presented, which should represent the time of pulse reflection at the water bottom in the full-waveform data. Figure 5 shows the options: threshold, peak maximum, centre of gravity, and zero crossing of the second derivative. With the exception of the centre of gravity detector (due to the influence of the water column), these standard methods are suited for the detection of both the water surface and bottom echo in ALB-based full-waveform data. Another widespread and more complex method of pulse detection in airborne laser scanning (ALS) data is Gaussian decomposition (Wagner et al. 2006;Reitberger et al. 2009;Mallet and Bretar 2009). The method is characterized by fitting several parameterized models to the measured receiver signal. However, Gaussian decomposition is only partially applicable to ALB data, since the signal shape can only be approximated to a limited extent by a Gaussian function due to the influence of the water column on the signal. For ALB applications, an approach called exponential decomposition was presented in Schwarz et al. (2017), which considers of the characteristics of bathymetric full-waveforms using exponential functions. This takes into account the influence of the water column in the fullwaveform. In Schwarz et al. (2019), the model is further specified by extracting the water surface, water column, and water bottom signals from the full-waveform. For our approach, a simple maximum detector is sufficient for detecting and extracting the water surface and bottom peaks as a first approximation.

Full-Waveform Analysis-Peak Significance
In general, a full-waveform may have several peaks or local maxima, some of which are due to interaction with objects but also to noise and erratic effects (Fig. 2). The peak detection (Sect. 3.2) returns all detected peaks unfiltered so that the relevant peaks must be determined by suitable methods. The goal is to detect the two most significant peaks, which are assumed to represent the water surface and water bottom. For this purpose, the peak description parameters isolation, prominence, and amplitude are used (Fig. 5). The parameters isolation and prominence have been adopted from Kirmse and de Ferranti (2017), where they are introduced for the detection of significant mountain peaks in elevation data. For the exponential decomposition in Schwarz et al. (2019), this approach is used to detect the most significant peak in ALB full-waveforms for initial value estimating.
The isolation describes the shortest distance between a local maximum and the signal part with the same signal strength. The prominence indicates the largest intensity difference in the isolation area, and the amplitude is here equal to the signal strength/intensity of the corresponding local maximum. Based on these parameters, the significance of each local maximum of the full-waveform can be estimated by multiplication of the three parameters. This approach has the advantage that local maxima evoked by signal noise or erratic effects can be filtered easily in most cases. In consequence, the second-most significant peak usually represents the water bottom.

Full-Waveform Analysis-Peak Width
The stacked full-waveform combines several bathymetric measurements with possibly different water depths or number of samples. Different water depths are due to roughness or slopes of the water bottom. To take these cases into account for later evaluation, the parameter peak width was introduced, which functions as a tolerance range around  6). This means that the search corridor is defined by the number of sample (water surface-bottom) and the peak width. The two parameters are determined using the stacked full-waveform.
It should be noted that by neglecting the terrain slope or underestimating the peak width, potential water bottom points could be lost. If the search corridor is too large, this can lead to ambiguities (more than one peak in the search corridor) and thus to incorrect detections of water bottom echoes in the individual full-waveforms. Water bodies with steep bottom topography require the consideration of additional information (Sect. 7).
The peak width is estimated by searching the local minimum between the bottom peak and the water surface peak that is closest to the bottom peak (Fig. 5). The number of samples between this minimum and the bottom peak's maximum is anticipated as the half peak width.
The peak width is used for filtering the full-waveform stacking results (Sect. 3.6) and for detection of water bottom echoes (Sect. 3.7).

Full-Waveform Stacking and Detection of Search Corridor
The joint processing of full-waveform data within a small spatial neighbourhood and consequently a similar laser pulse propagation path aims at a considerable reduction of signal noise and erratic backscatter effects within the water column. Herein, the difference in the forward and backward look of the Palmer scan pattern is not explicitly taken into account, but indirectly via the empirical determination of the processing unit size. In addition, the use of redundant information leads to an improved detection of weak water bottom echoes and consequently to a more reliable detection and extraction of deeper water bottom points. A faint water bottom echo candidate peak in an individual full-waveform is accepted if an according peak is present in the stacked full-waveform. This corresponds to a non-linear filtering, similar to the majority voting technique known from binary image processing. The basic requirement for combining full-waveforms is a correct alignment to each other with respect to the height coordinates (corresponding here to the z-coordinate; Fig. 6). This can be derived from the position of pulse emission, laser beam direction, and the laser pulse travel time. As this information was not available for our data, the maxima of the water surface echoes in the full-waveforms were used for vertical pulse alignment, with the assumption that the water surface approximately corresponds to a horizontal plane at the time of measurement. Subsequently, all measured full-waveforms signals of the respective processing unit are summed up resulting in a stacked full-waveform for every processing unit (Fig. 6). As mentioned above, the analysis of the stacked full-waveform is the basis for the determination of the search corridor for water bottom echoes.

Filtering of Full-Waveform Stacking Results
The full-waveform stacking is independently performed for each processing unit. The following ground echo detection works very reliably in stacked full-waveforms; but nevertheless, there are rare cases of incorrectly detected water bottom echoes. To handle these cases, the results of full-waveform stacking including the water bottom echo detection is checked iteratively in a filtering procedure using the information of neighbouring processing units. In every iteration step, all processing units, which are directly adjoined to processing units with reference water depth information or previously checked processing units, are searched and checked (Fig. 7). In detail, the approach averages out all checked or known number of samples (= water depths) as well as peak widths of adjoining processing units. This neighbouring information is then compared with the detected water depth or number of samples of the stacked full-waveform of the current processing unit. The use of the peak width enables coping with a certain slope or depth variation of the water bottom between the processing units. The following cases are possible, whereby the search corridor serves as a check corridor (based on the known information of the surrounding processing units): (1) the detected ground peak is within the search corridor; (2) the detected ground peak is outside the search corridor, consequently the most significant peak within the search corridor is the new ground peak; and (3) the detected ground peak is outside the search corridor and there are no peaks within the corridor range.
For cases (1) and (2), the processing units are considered as checked processing units, while for case (3), the processing units have no water depth information and are invalid for further iteration steps. This filtering mechanism is applied to all processing units classified as water.
To ensure reliable functioning of the procedure, the first iteration step requires processing units with reference water depths and the corresponding echoes in the full-waveform data. As reference values are not available in most cases, processing units with reliable water depth information have to be identified from the data set itself. Close to the shore, we know that the water depth is probably very shallow, but the determination of the water bottom echoes is not very reliable due to the superimposed water surface and water bottom echoes. Therefore, the use of the processing units near the shore poses a certain risk. In the middle of an unknown water body, water depth is likely to be too large for reliable water bottom echo detection. For this purpose, the parameters isolation, prominence, amplitude, and significance of the stacked full-waveform can help here as well. For reliable detection of water bottom echoes in full-waveforms, a clear separation (high isolation) and a clear distinction from the noise signal (high prominence) is advantageous. From this, it can also be deduced that a high significance indicates simple water bottom echo detection. To determine a suitable starting value of the water depth for filtering of the full-waveform stacking results, all stacked full-waveforms are used, and the parameters isolation, prominence as well as significance are analysed regarding their estimated water depth with the goal of finding either the water depth itself or at least processing units with the highest prominence and significance values. In detail, the mean of the prominence, isolation, and significance values are calculated as a function of their water depth. The starting water depth results from the maxima of the parameters. To consider that the stacked full-waveforms consist of a different number of measured full-waveforms, the isolation, prominence, and significance values are normalized regarding the water surface echo. The results of the start value estimation for the filtering algorithm are shown in Sect. 4.2.

Detection and Extraction of Bottom Point
Sections 3.5 and 3.6 focus on a reliable detection of the water bottom echo in stacked full-waveforms and subsequently the derivation of a search corridor. The final water bottom points are not derived from the stacked full-waveform itself as this leads to undesired smoothing effects (similar to a low-pass filtering) and a reduction of the bottom point resolution. Instead, the water bottom echoes are searched in the individual full-waveforms using the derived search corridor. In consequence, the estimated search corridor is used for the most likely water bottom echo in individual full-waveforms (Fig. 8a). Based on this information, the water bottom point of each individual full-waveform can be estimated.
First, the direction vector 1 (narrow laser beam) is calculated using the given points ls and ws (Fig. 8b). Then, the direction vector 2 is calculated by refraction of the direction vector 1 at the water surface according to Snell's law of refraction and the light transmission algorithm presented in Glassner (1989) (Eq. 1): with The water surface consists of measured water surface points and represents the optical interface between air and water. Due to a relatively high noise in the water surface points, a conversion into a regular grid was carried out to smooth the water surface. The laser beam is refracted at the smoothed water surface (air-water interface) using horizontal water (1) = n air ∕n water , n air ∕n water = refractive index air/water = water surface normal vector. Fig. 7 Iterative processing of a part of the study area. Brown = land grid cells. Yellow = currently checked water grid cells. Red = not yet checked water grid cells. Green = known or checked water grid cells.
Step 0 shows the unchecked study area.
Step 1 shows the water grid cells featuring reference or reliable starting water depth information.
The following iterative steps show the processing progress surface elements. The loss of accuracy compared to locally tilted water surface elements for low waves is negligible (see Sect. 5). Finally, the local maximum closest to the ground echo of the stacked full-waveform is used for calculation of the water bottom 3D-coordinates wb : where t 2 is the estimated travel time from the water surface to the water bottom (equal to the half sample time between the surface and ground echo) and c w is the light velocity in water.

Results
This section starts with the results of the determination of an appropriate grid cell size (used in this case as processing unit) and a suitable starting value for filtering the fullwaveform stacking. After that, the results obtained from processing data of the study areas are presented and evaluated to assess the quality and potential of the extended approach. For this purpose, the height deviation of the points of fullwaveform stacking processing (FWSP) to the points of ALB online processing (OWP) and echo sounder measurements (ES) are determined and used for accuracy analysis.
(2) wb = ws + t 2 ⋅ c w ⋅ 2 , Furthermore, the height deviations of FWSP points are compared against the total vertical uncertainty (TVU) of the Special Order released in the IHO Standards for Hydrographic Surveys (International Hydrographic Organization 2020): with The result analysis is carried out at first by visualization of the resulting point cloud to obtain first impressions regarding their plausibility and quality, followed by a quantitative comparison between FWSP points with existing water body bottom points (OWP and ES points). The quantitative evaluation contains the height deviation h , mean of height deviation h , sigma mean height deviation h , root mean square RMS, the sigma mean absolute deviation MAD(mean) , and the sigma median absolute deviation MAD(median) (Eqs. (4)-(9)) as well as reliability values that indicate the percentage of points not exceeding a given height difference: (4)  The height deviation and reliability values are presented for all points of the corresponding study areas, which allows a first indication of the quality of the FWSP data. However, these values vary notably with the water depth and are only representative for a part of the riverbed, which will become obvious in the following evaluation. In consequence, an investigation of the height deviations and reliability values as a function of their water depth is necessary. In addition, the results regarding the determination of the appropriate grid cell size and the start value for filtering the full-waveform stacking results, which are necessary for the best possible processing sequence, are briefly presented first.

Estimation of Grid Cell Size
To find an appropriate grid cell size for full-waveform stacking, different grid cell sizes from 1.5 m × 1.5 m to 3.0 m × 3.0 m were tested and compared regarding the number of extracted points, reliability and mean height deviation. The results are shown in Tables 2 and 3 whereby the boldprinted columns are absolute results and the other columns show the differences from the bold-printed columns.
A grid cell size of 2.0 m × 2.0 m (approx. 80-220 fullwaveforms per grid cell) delivers the best results for the investigated study areas Elbe 1 to 4 and Elbe 6 ( Table 2). For Elbe 5, a grid cell size of 2.5 m × 2.5 m (approx. 79 fullwaveforms per grid cell) was necessary due to a lower point density with 11 to 14 points/m 2 (Table 3). Investigations for this data set showed worse reliability results for smaller grid cell sizes and thus more incorrect estimated points. In contrast, using greater grid cell sizes, the number of FWSP points decrease considerably. The selected grid cell sizes thus represent a suitable compromise for good results.

Estimation of the Start Value for Filtering of Full-Waveform Stacking Results
To estimate a water depth for the filtering module (Sect. 3.6), the prominence, isolation, and significance of the water bottom candidates of all stacked full-waveforms were investigated as function of the water depth. As described in Sect. 3.6, high values for prominence, isolation, and consequently significance of an echo peak are good indicators for a reliable water bottom echo detection. Figure 9 shows the (a) isolation, (b) prominence, and (c) significance of large parts of the available data set (Elbe 1-Elbe 6), where the grid cells with a water depth of 0.75 up to 0.90 m have the highest values. Because an increasing water depth makes the separation of the water surface and water bottom echo more reliable, 0.9 m served as a starting value for the filtering algorithm. Figure 10a presents the resulting point cloud of the extended processing approach for the study area, where the points are colour-coded regarding their height coordinates. Compared to the data processed online, the results indicate an improved coverage of the riverbed (Fig. 10b) in larger depths. The other investigation areas feature similar results. The river cross sections in Fig. 11 confirm how well the FWSP points represent the actual bottom of the river. The following facts are already clearly visible in the visualization: Fig. 9 Results of evaluation of optimal start depth for full-waveform stacking filtering: a isolation, b prominence, and c significance as a function of their water depth. The black points are the single values, and the red line presents the mean trend. The vertical stripe pattern in the shallow water depths (0.3-0.9 m) is caused by the sample time interval (= resolution in laser beam direction of approximately of 6.5 cm) in the full-waveform data. At greater water depths, this effect is considerably weakened by the effect of the obliquely incident laser beam The visual impressions will be proven in the following by evaluating the FWSP data with regard to height deviation and reliability.

Analysis of Height Deviation
This section deals with the quantitative evaluation of the comparisons between FWSP points and reference points (OWP and ES points). The OWP points are mainly restricted to the shallower sections of the river, and the ES points mainly cover the deeper parts of the river due to vessel accessibility. The evaluation of both sets of comparison data is performed separately starting with the OWP points. Table 4 shows the results of the comparison with the OWP point cloud regardless of the water depth. According to the mean height deviations h , the FWSP points are slightly deeper (1.2-2.9 cm), but h is not significant in regard to h (significance level = 5% ). The RMS of 11-12 cm and the MAD(median) of 6.7-7.1 cm indicate small deviations across all investigation areas. In contrast, Fig. 13 shows a large variation of the height deviations between the OWP points and the FWSP points. The deviations in height seem to depend on the water depth. In particular, the large deviations in the transition zone from land to water are clearly visible. The FWSP points tend to be estimated too deep in this area. Therefore, Fig. 12 presents the statistical values (a) h , (b) RMS, (c) MAD(mean) and (d) MAD(median) as function of their water depth. The dependency between deviation measures and water depth is obvious. In detail, the points appear up to 0.3 m too deep in the shallower water section up to a water depth of 0.7 m. For water depths of approximately 0.7-1.8 m, the height deviations are small. In addition, the h values were mostly in the millimetre range so that the majority of the values shown in Fig. 12a were also significantly determined ( = 5%) and consequently more representative, as in Table 4. Furthermore, in Elbe 5, the deviations increase, starting from a water depth of 1.6 m. However, there are very few comparison points from OWP in these areas (blue dashed line in Fig. 12a-d), so this effect is difficult to prove. Table 5 presents the statistics of the comparison between FWSP and ES points regardless of the water depth. The FWSP points are determined approximately 1-5 cm too high. The h is not significant ( = 5% ) for the entire study area, similar to the comparison between FWSP and OWP points. The RMS values of 9.4-15 cm and the MAD(median) values of 8.7-9.9 cm are slightly larger than the comparison with the OWP points, but the FWSP points still have small height deviations. Figure 14 exhibits the height deviations of the FWSP points compared to the ES points, and again a large variation of height deviations is clearly visible. The evaluation as a function of the water depth is presented in Figure 15. Up to a water depth of 2.1 m, the deviation values are 10 cm or less. For a water depth of 2.2 m, the values are approximately 15 cm, but the number of detected FWSP points is already decreasing rapidly.

Reliability
For the evaluation of the reliability, the inlier rate estimation (Litman et al. 2015) was used, as the comparative data do not fulfil the requirements of ground truth data. The inlier rate indicates the portion of points whose height deviation does not exceed a given limit. As a limit orientation, the TVU value of the Special Order of the IHO for hydrographic surveying (0.25 m) was used. In addition, the evaluation was also performed for limit values ±10 cm around the IHO value.
Tables 6 and 7 show the number of extracted FWSP points that do not exceed a given maximum deviation from the comparison points (15 cm, 25 cm, 35 cm). Regarding the OWP data (Table 6) the FWSP points show a high inlier rate of 94.38% for a deviation of 25 cm or less. A   Fig. 11 Cross section of the Elbe River. Black = OWP points, Blue = water surface points, Red = ES points, Green = FWSP points maximum deviation of 15 cm is undercut by 88.80% of the FWSP points, and at a limit of 35 cm, the percentage of points increases to 97.51%. Similar results are obtained by comparison with the ES points (Table 7). For 15 cm, there was an inlier rate of 89.20 up to 99.34% of the points with a deviation of less than or equal 35 cm. As with the analysis of the height deviations, it is worth taking a look at the inlier rates in relation to water depth. Figure 16a-c shows the results of the reliability analysis with the OWP points, and Figure 16d-f shows those with the ES points. In general, a high inlier rate in areas with a water depth of 0.7 up to 2.1 m can be derived from all graphs, confirming the analysis of the height deviations. The problems in the transition area between land and water are also evident here. Furthermore, a strong decrease of the inlier rate can be observed from a water depth of 2.1 to 2.2 m, where approximately 11% of all extracted FWSP points were deeper than 2.1 m or approximately 2% deeper than 2.2 m. The considerable deterioration of the inlier rate as well as the already observed decrease in the number of extracted FWSP points suggest a maximum evaluable water depth of approximately 2.1 m, also confirming the height deviation analysis. It should be noted that this latter value only holds for the data set at hand, acquired in a rather turbid river.

Discussion
In the following, the results, observed problems, and additional information regarding the error budget are discussed to assess the results and to evaluate the potential of the fullwaveform stacking approach.

Comparison Results
Good results (RMS of 9-15 cm or MAD(median) of 6-10 cm) as well as high inlier rates (or reliability) were obtained for water depths from 0.7 cm to 2.1 m. At larger water depths, the height deviations and inlier rates became worse in our  pilot study data set with rather turbid river water, with the number of extracted points decreasing rapidly. The poor results in the range of water depths from 0.0 to 0.7 m can be caused by problems with the interpretation of superimposing echoes (scatterer cluster) from the water surface, water column and water bottom. This problem with superimposing echoes is also shown in Wagner et al. (2006) and Ullrich and Pfennigbauer (2011). In addition, the shape of the system full-waveform (shown in Schwarz et al. (2017)) leads to false detected water bottom echoes. In detail, the water bottom echo is superimposed either by the water surface or the volume backscatter. Consequently, the backscatter of the second peak of the emitted system full-waveform (with a low intensity) or another weak echo has a higher significance (Sect. 3.3) and is detected as water bottom echo. In this specific case, the ground points are mostly determined too deep in this area (Fig. 12a), attributed to a misallocation of the detected ground echo in the full-waveform. Commercial software packages, including those from the device manufacturers themselves, provide very good and plausible solutions in these cases so that the focus for the extended processing method remains in larger depths, where conventional methods show their limitations in finding correct water bottom echoes. Thus, a combination of the results from the standard OWP and from the full-waveform stacking will provide a good representation of the water bottom for water depths from 0.0 to 2.1 m in this case.

Error of Maximum Peak Detector
The use of the maximum peak detector for the detection of the water surface echo and the water bottom echo was simple to implement and suitable for first test purposes.  Furthermore, it already showed good accuracies compared to the echo sounder measurement. However, with this method, there is the risk that the real echo maximum is located between two samples so that for a detected maximum peak, the worst-case error must be assumed with the half sample width. Related to the specification in Sect. 2, the resulting height error for one detected peak is approx. 0.032 m in water and approx. 0.043 m in air. The exponential decomposition (presented in Schwarz et al. (2017Schwarz et al. ( , 2019) reflects the reality in a better way, and an increase in accuracy of the 3D coordinates of the ground points is expected.

Accuracy of Water Surface Points
Another important aspect that has been neglected in this article so far is the accuracy of the water surface point, which plays an elementary role in the consideration of the water bottom point accuracy. Equation 2 shows that the starting point for calculating the water bottom point is the water surface point. The water surface point in turn is indirectly obtained from the measured water surface points of OWP (Sect. 3.7), which are also subject to error influences due to the properties of the green laser. An empirical study Mandlburger et al. (2013) investigated the near water surface penetration of green laser signals and noticed a penetration depth into the water column of 10-25 cm. Furthermore, the authors stated that an underestimation of the water surface of 3 cm leads to a height deviation of 1 cm of the water bottom points. The derivation of the correct water surface could be achieved by statistical analysis and aggregation of neighbouring information.

Error Budget of the Water Surface Model
Moreover, the error budget caused by the wave pattern and the used water surface model has to be considered. With regard to the measured water surface points (e.g., in Fig. 11) the geometrical error budget of a water surface with smooth, rippled sea swell or small wavelets (wave height of 0.1-0.2 m) is considered. In Westfeld et al. (2017), the RMSE of lateral coordinate displacements ( XY ) and of the height coordinate displacement ( h ) is specified with 2.5-3.8% and 0.35-0.58% of the water depth respectively. Applied to the presented results, this means that for a water depth of 2.1 m, there is a lateral coordinate displacement of 0.053-0.08 m and a height error of 0.0074-0.012 m. To quantify the influence of the water surface model on the bottom point accuracy, further investigations were carried out with a bathymetry simulation tool with regard to the same wave types (Richter et al. 2018(Richter et al. , 2021. The results show that water surface modelling with horizontal water surface elements implemented in this study delivers small coordinate displacements at the water bottom with a lateral displacement of 1.0% of water depth and a depth displacement of 0.25% corresponding to 0.021 m and 0.0053 m at a water depth of 2.1 m respectively.

Conclusion
Airborne LiDAR bathymetry is an efficient area-wide measuring method for the determination of water bottom elevation models in shallow water areas. Due to water turbidity and the resulting limited penetration depth of the green laser beam into the water body, the full-waveform signal decreases with increasing water depth until a reliable detection of the water bottom in the measurement signals is no longer possible. To date, the maximum penetration depth is characterized by the reliable detectability of these echoes. In this article, an approach for the reliable detection of very weak water bottom echoes was presented with the goal to increase the penetration depth and therefore improve water bottom topography extraction. In contrast to the approaches focusing on isolated evaluation of full-waveforms, information from surrounding measurement data is used for this purpose, with the assumption that they show similar properties and characteristics. In the introduced non-linear full-waveform stacking techniques, closely neighbouring full-waveforms are combined to a stacked full-waveform, analysed, and the results are used for the determination of water bottom echoes in the individual full-waveforms to avoid smoothing effects. For the evaluation of the processing results, visualizations, an analysis of the height deviations, and inlier rates were used to show the potential of the presented approach. According to the shown results, a good representation of the riverbed can be assumed. In detail, the method shows small height deviations for water depths of 0.7-2.1 m ( RMS ≤ 14 cm , MAD(median) ≤ 8 cm ) and high inlier rates (94.38-97.90% with h ≤ 25 cm ). Based on the results of the evaluation, it can be stated that the penetration depth into the water body could be increased by 27.3% from approximately 1.65 m to approximately 2.1 m, which obviously leads to an improved coverage of the riverbed (approximately 204% related to the OWP covered area). It should be noted that the above quantitative statements are valid only for the present data set at hand, which was acquired in a rather turbid river.

Outlook and Future Work
In the work presented here, waveform stacking was applied to generate stacked full-waveforms, in which bottom echoes were detected by a technique analysing the peak parameter isolation, prominence and amplitude. A detected peak in an individual waveform was accepted (i.e., confirmed), if it was within a search corridor around the peak in the corresponding stacked waveform. This works well if there is only one local maximum (or peak) within this search corridor. To enhance the applicability to waterbodies with steep bottom topography, the technique will be extended by non-horizontal or non-planar search corridors obtained from neighbouring stacked full-waveforms.
In addition to the integration and evaluation of further methods to improve the shown approach, the applicability of the processing method for data sets from other ALB scanners and the potential for processing coastal data has to be investigated.
Regarding the adaptivity of the full-waveform stacking approach, a detailed investigation with clearly less neighbourhood information would be interesting because some bathymetry scanners acquire data with a considerably lower spatial resolution. In our own first experiments, data from the bathymetry scanner AHAB Chiroptera were processed with the full-waveform stacking approach, with an increase of the evaluable water depth up to 10%. The measurement resolution of the available data was less than 2 points/m 2 , requiring a considerably larger raster cell size. Another promising measuring system is the compact topo-bathymetric laser scanner Riegl VQ-840-G (RIEGL 2020). This measuring system is designed for use on an unmanned aerial vehicle, covers a square metre with approximately 20-50 points (Mandlburger et al. 2020) and achieves a maximum water penetration depth of twice the Secchi depth. The special feature of this laser scanner is the fact that for each emitted pulse, the complete full-waveform is recorded and stored, regardless of whether the online processing detects an object echo or not. Here, it is interesting to investigate whether object points can be detected and extracted from supposedly irrelevant measurement data.
Another aspect of future work will be the applicability of the methods to coastal waters, where lower turbidity and thus higher penetration depth can be expected. However, other conditions such as the impact of ocean waves to water bottom point estimation will represent new challenges. Therefore, the use of enhanced water surface modelling and consideration of wave-induced systematic errors are necessary. It is expected that the transfer of the presented processing methods to maritime ALB data will lead to an improvement in the acquisition of seabed topographies. etry is funded by the German Research Foundation (DFG). We also thank the German Federal Institute of Hydrology for the research assignment of a pilot study, which allowed intensive practical applications of the developed methods. In addition, we thank the Waterways and Shipping Office Dresden (WSA) for the provision of the survey data and the Airborne Hydro Mapping GmbH (AHM) company for the good cooperation in the preparation of the survey data. This research was supported by Deutsche Forschungsgemeinschaft (Grant 390023295).
Funding Open Access funding enabled and organized by Projekt DEAL.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.