Pre-Detection and Pre-Registration Averaging of Full Wave Signals in Airborne LiDAR Bathymetry

Abstract

A well-known technique to enhance the signal to noise ratio (SNR) of repetitive signals is to average them. The coherent parts of the signal add up constructively while the incoherent parts are averaged out. The prerequisite is that the signals are acquired under conditions of high repeatability, i.e., the signals must be sufficiently similar. In the present technical note, we describe an efficient method for maintaining signal similarity by ensuring spatial and temporal proximity of laser waveform signals obtained by a sensor operated from an airborne platform. The method makes use of a few auxiliary parameters such as laser pulse repetition rate, mirror rotation rate, platform altitude and flight speed. The method can be extended to be operated in real time.

1. Introduction

Light in the visible region has the ability to substantially penetrate into water bodies. A common wavelength of laser sources used for bathymetry is 532 nm which is visible as green light and is close to the optimum. Light detection and ranging (LiDAR) bathymetry measures the time a pulse of such green light needs to travel from the laser source to an underwater target and back to the instrument’s receiver. Furthermore, the power density is reduced through beam geometry and scattering. As the pulse propagates through air and water, it looses energy by interaction with the medium. As a result, the pulse amplitude at the receiver can become very weak, i.e., the signal-to-noise ratio (SNR) becomes very small, which can make reliable detection of the ground difficult.

One conceivable remedy to make the pulse stronger at the receiver would be to increase transmission power. Alternatively another method is to try to reduce the noise. On the one hand, it is conceivable to further improve the optical and electronic components, such as narrower optical filters and amplifiers that contribute very little noise. On the other hand, storing the received signal samples, a.k.a., Full-Waveform, enables SNR enhancement during post-processing.

Gutierrez et al. [1] mention the possibility of sensitivity enhancement by “averaging multiple waveforms over homogeneous surfaces” and call this stacking. Stacking generally refers to methods of averaging a number of measurements in order to enhance the coherent parts and attenuate the noise, thereby increasing the SNR. For a review of the history of stacking as used in seismic exploration, see Rashed [2]. Stacking methods already have been successfully used in LiDAR by Stilla et al. [3] who describe a method of how to automatically align neighboring waveforms in order to maximize coherency and detect weak pulses. They assume a plane geometric structure in the local neighborhood and sweep over a range of angles until a maximum is found. Mader et al. [4] use waveform stacking in LiDAR bathymetry. For this purpose, the origin and direction of the laser beam axes are first registered from the local sensor coordinate to the project coordinate system using the position information of the measuring platform. This makes it possible to use the water surface as a prominent geometric feature to align the single waveforms by shifting along their axes before averaging. For mainly horizontal ground geometries, they find that the recognizability of ground points is actually improved. It appears noteworthy that they do not directly use the stacked waveforms for estimation of the resulting target distances, but use the distance found from the stacked waveform to define a search range within which the final range is estimated from the single waveform.

Averaging, which is at the core of waveform stacking, has a long history in remote sensing. Motivated by a work of Lee et al. [5], George and Zamanakos [6] investigated so-called Comb Filters for enhancement of SNR in radar. Comb filters result from coherently adding up a periodic signal and can be implemented as a tapped delay line. This type of signal processing was called pre-detection averaging because the summation of the repeated pulses, taking into account the phase of the echo signal, takes place before signal detection, see, e.g., [7].

In LiDAR with the advent of digital signal processing, it became feasible to record the output signal of the optical detector. The coherent addition of several such signals before processing, i.e., detecting targets, is another closely related form of pre-detection averaging and has been used by Riegl [8] and Pfennigbauer et al. [9] for some time.

In this technical note, we will describe an efficient variant of pre-detection averaging that can be applied to full-wave laser data and discuss it specifically for use in bathymetry. In Section 2, we describe the datasets used for the assessment and explain the algorithms. In Section 3, we present the results and discuss them in Section 4. We conclude with Section 5.

2. Materials and Methods

2.1. Datasets

We used four datasets, two from the French Mediterranean coast covering the same area but flown with different parameters, a third from the Venetian lagoon and the fourth from the Venetian Porto di Lido-Canale di Nicolò (PdLCdN). Key parameters are listed in

Table 1. Summary of flight strip key parameters. Numbers are median values derived from trajectory data and registered point clouds; altitude is range as measured by online waveform processing (OWP) times cosine of nadir angle obtained from trajectory.

Site	Strips	Swath	Altitude	Speed	LPS	PRR	Density
		(m)	(m)	(m s−1)	(s−1)	(kHz)	(m −1)
Fr. 1 G1	9	47	76	2.8	101	199	467.1
Fr. 1 G2	4	46	75	19.9	101	199	69.4
France 2	12	46	76	2.4	97	51	149.6
Lagoon	25	100	161	56.8	100	51	2.9
PdLCdN	20	103	165	56.4	100	51	2.8

. The area at the French Coast was flown twice, with different pulse repetition rate (PRR) settings. The set with PRR = 200 kHz has been further grouped into two subsets because of the different ground speeds, the G1 set being closer to the shore. The datasets of Venice were recorded on 1 June 2024, and the French dataset on 19 December 2023. The latter was recorded from an uncrewed platform at a much lower flight speed. Swath and altitude data were calculated from online waveform processing (OWP) laser range results and off-nadir angles from trajectory data. Point density is a nominal value in the center of the swath, of either the forward or backward looking direction only. Swath width, however, only is a first-order approximation based on an average off-nadir angle, not addressing the ellipsoidal distortion of the scan pattern. The nominal value of mirror rotation speed was the same for all sets.

Although no data on water turbidity for the Venice Lagoon were collected during the survey, a value of 4 to 7 FNUs (FNU stands for Formazin Nephelometric Unit, see [10,11]) at the day of data acquisition can be found from satellite observations using the European Union-Copernicus Marine Service [12] which is given, however, outside the lagoon. At the day of acquisition of the French datasets, no value is available from European Union-Copernicus Marine Service [12], but at 25 December 2023 a reading of 0.6 FNUs and at 25 December 2023 a reading of 0.3 FNUs is available from the same service.

2.2. Sensor

All Datasets were acquired by sensors of type VQ-840-G from company Riegl (Horn, Austria) [13]. This sensor type is a compact topo-bathymetric laser scanner for airborne platforms with OWP and full waveform recording for full waveform analysis (FWA) in post-processing. The recording is unconditional, i.e., no trigger condition depending on presence of a target is required to initiate the acquisition of data samples. This is an essential prerequisite for using the method described in this technical note, as the signal before averaging may be too small to trigger the recording.

The laser wavelength of the instrument is in the green at 532 nm with a selectable beam divergence and receiver field of view. The scanner is equipped with a rotating mirror that is arranged in such a way that the nadir angle of the resulting laser beam has little variation over a scan line. The effective scan pattern is nearly elliptic with an off-nadir angle of $\pm {20}^{\circ }$ perpendicular to flight direction and $\pm {14}^{\circ }$ in flight direction. The pulse repetition rate is selectable between 50 kHz and 200 kHz, while the angular rotation of the mirror is between 10 and 100 lines per second.

The sensor was also equipped with internal inertial measurement unit and global navigation satellite system (IMU/GNSS) unit.

2.3. Pre-Detection Averaging

If N received signals

$$\begin{array}{c}\hfill \begin{array}{c}\hfill {\mathit{r}}_i\left(t\right)=p\left(t\right)+{\mathit{n}}_i\left(t\right),\end{array}\end{array}$$

(1)

i the signal number where $0<i\le N$ and $t=0$ the time of emission are modeled as disturbed by additive uncorrelated Gaussian noise with with zero mean and variance ${\sigma }^2\left(t\right)$,

$$\begin{array}{c}\hfill \begin{array}{cccc}\hfill \mathrm{E}\left\{{\mathit{n}}_i\left(t\right)\right\}& =0\mathrm{and}\hfill & \hfill \mathrm{E}\left\{{\mathit{n}}_i\left(t\right){\mathit{n}}_j\left(t\right)\right\}& =\left\{\begin{array}{cc}0\hfill & i\ne j\hfill \\ {\sigma }^2\left(t\right)\hfill & i=j\hfill \end{array}\right.,\hfill \end{array}\end{array}$$

(2)

the pre-detection averaging estimator $\hat{\mathit{p}}\left(t\right)$ of $p\left(t\right)$ is defined by

$$\begin{array}{c}\hfill \begin{array}{c}\hfill \hat{\mathit{p}}\left(t\right)={\displaystyle \frac{1}{N}}\sum _{i=1}^N{\mathit{r}}_i\left(t\right)\end{array}\end{array}$$

(3)

which has a smaller variance than estimating $p\left(t\right)$ by a single instance of ${\mathit{r}}_i\left(t\right)$. The second-order statistical properties of this estimator are

$$\begin{array}{c}\hfill \begin{array}{c}\hfill \mathrm{E}\left\{\hat{\mathit{p}}\left(t\right)\right\}=p\left(t\right)\end{array}\end{array}$$

(4)

for the expected value and

$$\begin{array}{c}\hfill \begin{array}{c}\hfill \mathrm{E}{\{\left(\hat{\mathit{p}}\left(t\right)\right)-p\left(t\right))}^2\}={\displaystyle \frac{1}{N}}{\sigma }^2\left(t\right),\end{array}\end{array}$$

(5)

for the variance, which is reduced by the factor N, the number of averaged signals. The SNR, a measure of statistical significance, is given by

$$\begin{array}{c}\hfill \begin{array}{c}\hfill \mathrm{SNR}\left(t\right)={\displaystyle \frac{\left|\mathrm{E}\right\{\hat{\mathit{p}}\left(t\right)\left\}\right|}{\sqrt{\mathrm{E}\left\{{\left[\hat{\mathit{p}}\left(t\right)-p\left(t\right)\right]}^2\right\}}}}=\sqrt{N}{\displaystyle \frac{\left|p\right(t\left)\right|}{\sigma \left(t\right)}}\end{array}\end{array}$$

(6)

and can be seen to increase proportional to the square root of the number of averaged signals.

2.4. Depth Performance

Maximum depth measurable with LiDAR depends on water turbidity. The Secchi depth $Z_s$ as standardized in ISO 7027-2 [11] is defined as “the depth at which a matt, circular disc no longer is visible”. The method, although not providing an exact measure of water clarity, since influenced by environmental factors, e.g., sun’s glare and operators eyesight capabilities, is a commonly used gauge for illustrative indication of performance in laser bathymetry.

The signal power at the receiver $P_r$, reflected from the water bottom at depth z, follows an exponential law

$$\begin{array}{c}\hfill \begin{array}{c}\hfill P_r=P_t\kappa e^{-2\gamma z}\end{array}\end{array}$$

(7)

with transmitted peak power $P_t$, system factor $\kappa$ and attenuation $\gamma$. The attenuation factor $\gamma$ and the Secchi depth $Z_s$ are related by

$$\begin{array}{c}\hfill \begin{array}{c}\hfill \gamma Z_s=\eta \end{array}\end{array}$$

(8)

where $\eta$ is reported in the range 1.1 to 1.7 by Guenther [14], 1.7 to 2.3 in ISO 7027-2 [11] or 1.48 in a recent publication by Lee et al. [15].

Assuming the minimum SNR required to detect a target at depth z is

$$\begin{array}{c}\hfill \begin{array}{c}\hfill \mathrm{SNR}={\displaystyle \frac{P_t\kappa e^{-2\gamma z}}{\sigma }},\end{array}\end{array}$$

(9)

to detect a target at $z+\Delta z$ the following equation can be established by using Equation (6)

$$\begin{array}{c}\hfill \begin{array}{c}\hfill {\displaystyle \frac{P_t\kappa e^{-2\gamma z}}{\sigma }}=\sqrt{N}{\displaystyle \frac{P_t\kappa e^{-2\gamma (z+\Delta z)}}{\sigma }}.\end{array}\end{array}$$

(10)

Canceling common factors on both sides the following expression can be obtained indicating the expected averaging count N required to increase the resolvable depth by $\Delta z$

$$\begin{array}{c}\hfill \begin{array}{c}\hfill N=e^{4\eta {\displaystyle \frac{\Delta z}{Z_s}}}.\end{array}\end{array}$$

(11)

2.5. Coherence Enhancement

Equation (1) is applicable in a strict sense only as a model of a laser echo signal in a non scanning situation. In a scanning situation every ${\mathit{r}}_i\left(t\right)$ will be a received signal from a different position and direction of the laser beam. Consequently, Equation (1) must be changed to

$$\begin{array}{c}\hfill \begin{array}{c}\hfill {\mathit{r}}_i\left(t\right)=p\left(t\right)+\Delta p_i\left(t\right)+{\mathit{n}}_i\left(t\right)\end{array}\end{array}$$

(12)

where it is assumed that the $\Delta p_i\left(t\right)$ are deterministic differences caused by the difference in target geometries due to scanning. $\Delta p_i\left(t\right)$ is zero for the one measurement treated as reference. In consequence, estimator (3) becomes biased

$$\begin{array}{c}\hfill \begin{array}{c}\hfill \mathrm{E}\left\{\hat{\mathit{p}}\left(t\right)\right\}=p\left(t\right)+b\left(t\right),\end{array}\end{array}$$

(13)

with bias

$$\begin{array}{c}\hfill \begin{array}{c}\hfill b\left(t\right)={\displaystyle \frac{1}{N}}\sum _{i=1}^{N-1}\Delta p_i\left(t\right).\end{array}\end{array}$$

(14)

The variance of the estimator is biased too, viz

$$\begin{array}{c}\hfill \begin{array}{c}\hfill \mathrm{E}\left\{{\left[{\hat{\mathit{p}}}_i\left(t\right)-p\left(t\right)\right]}^2\right\}=b^2\left(t\right)+{\displaystyle \frac{1}{N}}{\sigma }^2\left(t\right).\end{array}\end{array}$$

(15)

Enhancing coherence for a scanning instrument operated on a moving platform means reducing the bias by making the laser beam origins as close and their directions as parallel as possible. The circular scan pattern of Figure 1 is typical for bathymetric scanners for it allows to approach an almost constant nadir angle of the beam. Incremental movement $\Delta x$ of the platform is superimposed with an incremental rotation $\Delta y$ of the laser deflection mechanism. A number of adjacent shots, although close in time, can be spread out spatially quite substantially, as is indicated by the red circles in Figure 1. For the same number of laser shots, the maximum spatial extent is reduced by the grouping of shots not only in time but also in space. The green circles in Figure 1 illustrate the idea.

The nominal spot size of the laser when hitting the ground is the diameter of the green or red circles. The virtual spot size is then defined as the largest distance within a group of laser shots that hit the ground and are designated for averaging. Obviously, the green configuration, approximating a square, results in a smaller virtual spot size. The optimal number of lines $n_x$ and shots $n_y$ for a given averaging number is determined by searching for factors $n_x$ and $n_y$ such that $n_x\Delta x\approx n_y\Delta y$ while keeping the product $n_x{n}_y$ close to the given averaging number. For post-processing, the required values $\Delta x$ and $\Delta y$ can be obtained from the flight trajectory and measurement data by

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill \Delta y& =r\sin \left(\theta \right)\Delta \varphi =r\sin \left(\theta \right){\displaystyle \frac{\omega }{\mathrm{PRR}}}\hfill \\ \hfill \Delta x& =v{\displaystyle \frac{2\pi }{\omega }},\hfill \end{array}\end{array}$$

(16)

with r the average distance to the water surface, $\theta$ the nadir angle, $\Delta \varphi$ the angular shot increment, v the speed of the platform and $\omega$ the angular rotation speed of the scan mechanism. $\Delta \varphi$ is related to $\omega$ by the PRR. As these data are also available during the flight, i.e., prior to point cloud registration, the method described in this note can in principle also be carried out in real time. To be precise, Equation (16) is only exact for a circular scan pattern and even then only in the center of the swath. Nevertheless, we use it as a first-order approximation. We emphasize that the above reasoning is in the local sensor coordinate system, i.e., before registering the scan pattern to the project coordinate system.

In general, the times of flight of two echo pulses that occur at different beam angles on an inclined plane are unequal. While in theory the beam angles with respect to a horizontal water surface or a horizontal water bottom should be the same for different shots in a quasi-circular pattern as in Figure 1, this is not the case in practice. Pitch and roll angles of the platform as well as non-horizontal ground cause deviations from the ideal situation, and thus different times of flight.

Figure 2 shows the general case. Two echos from laser beams A and B, enclosing angle $\alpha$, are reflected from a slanted plane w with the angle $\beta$ w.r.t. beam direction O–A. Assuming that the velocity in direction of flight is small enough so that the platform does not move much between successive laser shots, a fixed point O is used for both beams. Neglecting pulse spreading due to beam divergence and incidence angles on the plane, the two echoes will be of similar shape. If both angles are known the pulse from beam B can be stretched by scaling

$$\begin{array}{c}\hfill \begin{array}{c}\hfill p_B^{\prime }\left(t\right)=p_B\left(t{\displaystyle \frac{\sin (\beta -\alpha )}{\sin \left(\beta \right)}}\right),\end{array}\end{array}$$

(17)

aligning the pulse location with pulse $p_A\left(t\right)$ at the cost of shape distortion. Differently to Stilla et al. [3], who perform a search for the optimum angle $\beta$ which maximizes coherence, we rely on the angle $\alpha$ being small. The same grouping criterion, which guarantees spatial proximity of the laser shots and was explained above, also leads to approximately parallel beams and thus small angles $\alpha$.

The horizon of the measuring platform is known, as are the angles $\alpha$; thus, it is possible to correct the pulse scales for special planes such as the water surface or the ground. We omitted such correction, however, based on the following consideration. A first-order Taylor approximation of Equation (17) for small $\alpha$

$$\begin{array}{c}\hfill \begin{array}{c}\hfill p_B^{\prime }\left(t\right)=p_B\left(t(1-\alpha cot\left(\beta \right))\right)\end{array}\end{array}$$

(18)

reveals that the relative stretching error is proportional to $\alpha cot\left(\beta \right)$. For $\beta ={90}^{\circ }$ the error vanishes and for $\beta ={45}^{\circ }\mathrm{to}{135}^{\circ }$ the $|cot(\beta \left)\right|$ is smaller than $1.0$ meaning for angles $\alpha <\pm 1^{\circ }$ the relative error will be 1.7% at most.

2.6. Sliding Buffers

Processing starts by splitting the data of each flight strip from the quasi-circular scan pattern in Figure 1 into two partial strips, one strip in the forward direction and one strip in the backward direction. The splitting defines scan lines: A scan line comprises the laser shots from one edge of the swath to the other. Each of the two sub strips is then processed separately so that the beam directions of neighboring laser shots are approximately the same, ensuring the small angle $\alpha$ condition stated in the previous section.

With the number of lines $n_x$ and the number of shots $n_y$ used for an averaged waveform, the total number of scan lines held in memory at any one time is always equal to $n_x$.

The data of the full waveforms are loaded into a buffer structure for processing which can be addressed via the line number, the shot number and the sample index as sketched in Figure 3. A batch of $n_x$ lines is loaded for a processing cycle. After processing a batch of lines, the buffer space of the oldest line can be reused when loading the next line, as the data of the oldest are no longer needed. A simple copy operation is not sufficient, however, for the sample values as they must be brought into a common grid because the relative sampling phases of the waveforms will be different in general. We use a simple linear interpolation scheme to resample the waveforms into a common grid of sampling instances. Since the required operation also can be seen as a shift of one waveform with respect to the other an implementation as a fractional delay filter [16] is possible under the assumption that the signals are bandlimited.

All sample values indicated by the dashed rectangle in Figure 3 contribute to an averaged sample $S_{\mathtt{avg}}$. We do not, however, perform a plain averaging operation, but also do an outlier removal by means of a simple percentile filter which effectively removes spurious samples that can occur as a result of the very sensitive sensor. This filter can be switched off at the user’s discretion.

A cycle of the averaging process finally runs in such a way that all lines of the outer cuboid are loaded in which the gray sub-cube then slides along the shot axis and at each position of the gray sub-cuboid the innermost rectangle slides along the sample axis where each position results in a single averaged sample $S_{\mathtt{avg}}$. The process continues with the next cycle by loading the data of the next line until all lines have been processed.

The total number of averaged waveforms will thus be close to the number of raw waveforms. A small difference results because at both ends of the swaths from the forward and backward looking sub-strips no averaging is possible.

2.7. Bathymetry Map and Image Preparation

The averaged waveforms were analyzed using the Surface Volume Bottom (SVB) algorithm by Schwarz et al. [17] and converted into a globally registered point cloud. After classification of the water surface points, a water surface model was created and used to apply refraction correction to the underwater points. The classification of the bathymetric points in the turbid Venetian datasets was done manually. The classification of the bathymetric points from the French datasets was done by using a depth threshold to cut off the water surface followed by some manual corrections.

RiProcess (Vers. 1.9.5)/RiUNITE (pre. v.1.0.8) from the company Riegl [18] was used for averaging, SVB analysis and refraction correction, while the free software tool QGIS (v.3.36.3) [19] was used for map preparation and 3D visualization.

3. Results

The nominal averaging count for post-processing of all sites was set to 100. The criterion for obtaining a shape as square as possible led to virtual spot diameters of ≈1.2–2.2 m for the French datasets and ≈8.8–9.2 m for the Italian datasets. The virtual diameter is the sum of of the intrinsic spot size and the diagonal’s length of the averaged area. The smaller virtual spot size of the French datasets is due to a a higher point density resulting from the slower airspeed of the uncrewed platform. A summary of parameters is given in

Table 2. Overview of site specific processing parameters. Numbers are median values, virtual spot diameter is the sum of beam spot diamter and the diagonal length of the averaged area. Aspect ratio $\Delta y/\Delta x$ is calculated from sensor parameters, i.e., mirror rotation lines per second (LPS) and PRR. N and averaged size follow from evaluation of the squareness criterion for selection of the lines and shots subject to averaging.

Site	Virt. Spot Dia.	$\Delta \mathit{y}/\Delta \mathit{x}$	N	Averaged Area
	(m)			(m × m)
Fr. 1 G1	1.2	2.7	98	0.5 × 0.4
Fr. 1 G2	2.2	0.4	96	1.2 × 1.2
France 2	1.6	11.3	93	0.8 × 0.8
Lagoon	8.8	1.1	108	6.3 × 6.1
PdLCdN	9.2	1.1	110	6.2 × 6.3

. N is given by the product $n_x\times n_y$ and the averaged area is $n_x\Delta x\times n_y\Delta y$. It can be seen that the averaged area always is almost square despite different point densities in flight and across directions.

The effectiveness of pre-detection averaging is demonstrated in Figure 4 and Figure 5. Figure 4 showing two randomly selected pulses from the neighborhood of a shot, in thin gray, and the average of ≈100 surrounding pulses, in solid black. The averaged pulse is of smoother shape with a much lower number of local extrema. The numbers of averaged pulses, however, is only approximate for two reasons: first because the actual number is subject to the squareness criterion as described in Section 2.5, second because outliers are removed prior to averaging.

For Figure 5, two profiles from France 2 dataset were taken at the same crosscut location. The upper profile is without averaging, while the lower profile was created using an effective averaging number of $N=93$. The maximum resolvable depth in the not averaged case is 19 m and 24 m with averaging. The improvement due to averaging thus is 5 m.

Although turbidity in the Venice area was high, underwater points could be determined from all three areas surveyed. Figure 6 shows a map of Porto di Lido-Canale di Nicolò superimposed with bathymetric points mapped to color. The maximum penetration depth of the laser below the water surface was about 6 m on the day of recording. The underwater structures of the MOSE (Modulo Sperimentale Elettromeccanicao), the Venetian flood protection system, are not visible in the bathymetry because when it is not in operation its parts lie below the penetration depth of the laser.

Maximum penetration depth within the Venetian Lagoon was similar to that of the PdLCdN dataset; the color map, however, was set to cover 3 m only to emphasize the high capability of structure resolution in the shallow water. Figure 7 exhibits a map of underwater bathymetry of the surveyed area. Some voids visible in the north–south direction are due to insufficient overlapping of flight strips. The map consists of bathymetric underwater points only and shows noisy gaps in places where land cover would be expected, as the model of the SVB algorithm is only suitable for underwater analysis. No maps for comparison purposes could be created without prior averaging of the datasets, because classification of bathymetric points was severely hindered by the large amount of spurious points in that case.

Turbidity was very low at the French Mediterranean Coast and the data showed that a depth penetration of 24 m could be attained. Figure 8 shows a detail view of the water bottom in perspective projection, taken from the France 1 G1 dataset, highlighting that it is possible to obtain structures in high resolution despite a high averaging number provided flight conditions are appropriate. The rocks in the centre axis of the picture are of diameter approximately 6 m. The visual appearance of the rocks was further enhanced by applying a shadowing effect. Maximum depth could not be achieved in the France 1 dataset, because of a storage limitation that applies for long waveform recordings taken with a setting of 200 kHz for the PRR.

4. Discussion

The effect of summation of a group of closely spaced laser beams with small footprint can be understood as being equivalent to a single pulse of higher energy but larger footprint. Consequently, it can be expected that the measurement properties of the presented method with respect to environmental conditions in terms of wind speed or water currents are comparable to a system with correspondingly higher laser power and larger beam divergence. It is known, e.g., that proper selection of beam divergence, and thus footprint size, can be used to compensate for the influence of surface waves using the averaging effect over the footprint size, see, e.g., [14]. Although the purpose of averaging is different in that case, i.e., not to increase the SNR, the principle is the same, which in our opinion is a further justification for the term virtual footprint caused by the group of laser beams subjected to averaging.

As indicated in Equation (6) the SNR should increase in proportion to the square root of the number of waveforms averaged. Instead of showing the quantitative improvement of the SNR in a direct manner, we use the reasoning from Section 2.4 to show the improvement indirectly. If both the increase of SNR by the square root of averaged waveforms and the exponential law of power decay are true, then the improvement in depth penetration, following from Equation (11) must be logarithmic in the number of averaged waveforms:

$$\begin{array}{c}\hfill \begin{array}{cccc}\hfill \Delta z& ={\displaystyle \frac{lnN}{4\eta }}{Z}_s\mathrm{i}.\mathrm{e}.,\hfill & \hfill \hspace{1em}\hspace{1em}\hspace{1em}5 \mathrm{m}& \approx {\displaystyle \frac{ln93}{4\times 2.15}}\times 9.5 \mathrm{m}.\hfill \end{array}\end{array}$$

(19)

From experience [20], we learned that typically two times the Secchi depth can be reached with the particular sensor type with SVB waveform processing. In an attempt to compare the theoretical relations with the empirical findings without having available on-site measurements of Secchi depth, we use the top profile in Figure 5, where a maximum depth penetration of 19 m is achieved, to estimate $Z_s$ to be 9.5 m. Furthermore, a value of 2.15 for $\eta$ is taken in accordance with ISO 7027-2 [11]. With these values, relation (19) is met quite closely.

One of the adverse effects of averaging is, as must be expected, a reduction in resolution of the target scene due to the low pass behavior of averaging. This effect, however, depends on the sampling parameters of the un-averaged dataset as can bee seen in Figure 8. The strong oversampling in the France 1 G1 dataset still allows a high averaging rate without compromising resolution. But even if resolution is reduced, a certain degree of averaging will be highly beneficial because robustness against the influence of surface waves is increased, thus improving derivation of the water surface. Theoretically, a low pass filtered dataset should be able to be represented with less data without loosing information, as stated by the Nyquist-Shannon sampling theorem. We took, however, the approach of not reducing the amount of data, since often the higher resolution dataset is required for further processing anyways. On the other hand, since the storage and full waveform processing requirements can be demanding the potential for optional reduction may be interesting for some projects.

Even if a reduction in resolution is undesirable, a small averaging degree can still significantly reduce spurious targets, which are more likely to occur in turbid water. While these targets may be of interest for determining water properties in general, they make it difficult to efficiently classify the ground. Applying the proposed method during post-processing allows it to be adapted to different problem definitions.

Although the virtual footprint size of the Italian datasets is already half the size of satellite based systems, e.g., for ICESat-2 the footprint size is 17 m [21], the airborne system benefits from tunability of the size and can thus be adapted during post-processing. Similar to the approach of Mader et al. [4] one can also combine bathymetry with high averaging and one without averaging to increase the resolution by considering targets close to the smoothed bathymetry only.

Another possibility is to increase depth performance without reducing the resolution by scanning faster and flying slower yielding a higher point density an approach particularly interesting for use with uncrewed platforms. Figure 8 is an example of such a use case.

All test datasets show that the approximate Equation (16) suffices for selecting a suitable group of waveforms for averaging. A more thorough analysis would start at a more accurate description of the laser pattern such as the Palmer Scan, i.e., an elliptical pattern, see [22], or at an even more advanced analysis as conducted by Li [23] taking care of the deviations from the elliptic shape.

Another question that requires consideration is the determination of the mean distance r used in Equation (16). We use the median value that can be obtained from online range measurements, i.e., before the pre-detection step. This means that although these values represent plausible distances above a water surface, they do not necessarily reflect the distance to a water surface in heterogeneous areas with mixed water and land cover. For the Lagoon dataset, a discrepancy between the altitude above the water surface as stated in Table 1 and the surface as obtained from the registered pointcloud, as in Figure 7, was observed; however, the effect was not very large. It is a simple extension to use a nominal, manually specified, altitude instead, but at the expense of non-automatic operation.

We described the method for use in post-processing, but an extension that allows averaging during the flight in real time is conceivable. From the parameters needed in Equation (16), everything is available in real time from the IMU/GNSS unit.

Comparing the presented method to an approach where the waveforms are registered into a geographical frame before averaging, we can identify the following: Proper registration requires a sufficient number of points that can be derived before averaging is applied. The fact that two points are close to each other in geo-referenced space does not mean that the laser directions associated with them are parallel, a prerequisite for proper averaging.

It should not go unmentioned that the method allows a high degree of parallelism in processing and at the same time has a relatively small, constant memory requirement. This makes it possible to assign each of the horizontal slices that contribute to the $S_{\mathtt{avg}}$ point in Figure 3 to a separate CPU for simultaneous processing.

While increasing depth performance is one of the goals of averaging, within the scope of this technical note we could also demonstrate the effect of better structure recognition in shallow water. A higher SNR reduces the amount of spurious points resulting in a cleaner pointcloud which greatly assists point classification as demonstrated by the results of Figure 6 and Figure 7 which would not have been possible with the un-averaged datasets.

The noise model used in Section 2.3 is the ubiquitous Gaussian model under assumption of uncorrelatedness of samples. For a more thorough analysis, however, the model should be refined. A distinction would need to be made between sources of interference in the optical and electrical regimes. In the optical case, under high turbidity, the assumption that the repeatability of the pulse shape under stationary conditions is sufficiently high becomes questionable. In the electrical regime, the effect of the detector, a magnitude detecter, is to change the Gaussian noise into a Rician distribution. Since the signals are also amplified, another source of noise that depends on the magnitude of the signal, i.e., shot noise, enters into the refined analysis.

5. Conclusions

In this technical note, we describe a variant of waveform stacking, a method for increasing the SNR called pre-detection averaging, which is particularly useful as a pre-processing step of full waveform analysis to derive bathymetry. The higher SNR makes possible higher depth resolution and better point classification due to the smaller amount of spurious points, however at the cost of spatial resolution if the original dataset before averaging is not dense enough. The presented method is fast, requires only a small and constant amount of memory, can be applied to the unregistered dataset and can optionally be performed in real time during flight. The single parameter that needs to be set manually is the degree of averaging. The method can be understood as a means of processing the data with a virtual footprint size that is larger than the native footprint size. When carried out during post-processing, the optimum averaging size can be tuned for best results without the need for new data acquisition.