WATER TURBIDITY ESTIMATION FROM LIDAR BATHYMETRY DATA BY FULL-WAVEFORM ANALYSIS-COMPARISON OF TWO APPROACHES

Airborne LiDAR bathymetry is an efficient technique for surveying the bottom of shallow waters. In addition, the measurement data contain valuable information about the local turbidity conditions in the water body. The extraction of this information requires appropriate evaluation methods examining the decay of the recorded waveform signal. Existing approaches are based on several assumptions concerning the influence of the ALB system on the waveform signal, the extraction of the volume backscatter, and the directional independence of turbidity. The paper presents a novel approach that overcomes the existing limitations using two alternative turbidity estimation methods as well as different variants of further processed full-waveform data. For validation purposes, the approach was applied to a data set of a shallow inland water. The results of the quantitative evaluation show, which method and which data basis is best suited for the derivation of area wide water turbidity information.


INTRODUCTION
Airborne LiDAR bathymetry (ALB) allows the efficient surveying of shallow water bodies (Mandlburger, 2020). Laser bathymetry beam propagation is characterized by scattering and absorption effects in the water column, resulting in a continuous decrease of the received signal intensity (Churnside, 2014). The degree of intensity decrease mainly depends on the turbidity of the water body (Walker and McLean, 1999). Inverting this relation, statements on the local level of turbidity in the water body may be obtained by analyzing the decay of the recorded waveform signal.
Area-wide information on water turbidity provides an essential input for a wide range of ecological studies (e.g. Devlin et al. (2008); Gippel (1995)). Common techniques for turbidity measurement can be divided into in-situ methods (e.g. measuring of visibility depth with a Secchi-disk (Secchi, 1864), measuring of down welling irradiance with submersible quantum sensors, measuring of a beam attenuation coefficient with transmissometers) and laboratory methods (analyzing fixed volume samples with absorptiometers or nephelometers). The obvious disadvantages of all these techniques are the large effort, the high weather dependence, and the point-wise character of the measurements. Phillips et al. (1984) and Billard et al. (1986) present an approach for the area-wide derivation of water turbidity from fullwaveform (fwf) data of a large footprint profiling deep water ALB system, designed to capture ocean waters. They use the amplitude and the decay of the backscattered signal to determine the absorption and scattering coefficients of sea water. The * Corresponding author determination of the scattering coefficient requires a system calibration with in-situ reference data.
In Richter et al. (2017), we presented a first approach on the determination of area-wide turbidity information from smallfootprint full-waveform airborne LiDAR bathymetry data. The method is based on the analysis of the recorded waveform signals by fitting an exponential function. Herein the exponential decay coefficient depicts an integral measure describing turbidity. The results are very promising and show a high correlation with turbidity measurements obtained by conventional techniques. However, the approach is limited by the following aspects: (a) Using the signal part after the maximum for exponential function approximation may lead to a bias of the turbidity parameters due to dominant water bottom echoes.
(b) Deriving turbidity parameters from received waveforms leads to a distortion of the turbidity parameters due to the influence of the system waveform.
(c) The turbidity parameters obtained from single fullwaveforms signals may be strongly affected by signal noise.
(d) The signal propagates obliquely through the water column, i.e. in contrast to common techniques, turbidity is determined in the direction of the laser pulse propagation and not vertically downwards.
In this contribution, we want to overcome these limitations and therefore present a refined approach that uses only the volume backscatter segment of the full-waveform signal for turbidity derivation (a). Additionally, the surface-volume-bottom (SVB) algorithm presented in Schwarz et al. (2019) has been implemented, which was developed to reconstruct the differential backscatter cross section (dBCS) from bathymetric waveforms. The dBCS allows for the determination of turbidity information free from influences of the system waveform (b). To handle problems with signal noise and to enable a vertical derivation of turbidity parameters (c and d), we combine individual full-waveform signals to stacked full-waveforms  as well as ortho full-waveforms (Pan et al., 2016) resulting from a voxel-based data representation.
Our study compares the performance of the two approaches exponential function approximation (Section 2.1) and SVB algorithm (Section 2.2) to estimate the turbidity. For this purpose, we use individual full-waveforms, stacked full-waveforms, and ortho full-waveforms (Section 3) of a shallow inland water data set. In Section 4.1 the results are presented and discussed. Since no ground truth data are available, the comparison is based on a quantitative evaluation of the derived turbidity parameters. The evaluation is described in Section 4.2. The contribution ends with the conclusion in Section 5.

WATER TURBIDITY ESTIMATION METHODS
Bathymetric waveforms result from the interaction of the laser signal with the water surface, the water column and the bottom of the water body (Abdallah et al., 2012;Pfeifer et al., 2016). Figure 1 shows a schematic representation of the signal components. Analyzing the signal decay and the amplitude of the water column segment of the waveform allows for statements on turbidity. Essentially, the volume exponential decay coefficient k approximates the volume absorption coefficient a for large receiver apertures (Guenther et al., 2000). Moreover, it is almost identical with the diffuse attenuation coefficient (Gordon, 1982). The amplitude of the volume backscatter is indicative of the scattering coefficient b. Since the exact determination of b requires a calibration of the ALB system (Muirhead and Cracknell, 1986), it is not the subject of this work. Sections 2.1 and 2.2 explain the methods used to derive the volume exponential decay coefficient k.
water surface echo water bottom echo signal intensity sample / time volume backscatter Figure 1. Schematic full-waveform representing the three signal components according to Guenther et al. (2000).

Exponential function approximation
First, the volume backscatter must be extracted from the fullwaveform. For this purpose, the volume backscatter is defined as the signal part after the maximum of the water surface echo to the local minimum before the water bottom echo. The detection of the water surface echo and the water bottom echo is carried out by the full-waveform processing approach presented in Mader et al. (2019Mader et al. ( , 2021. Subsequently, an exponential function is fitted into the volume backscatter signal, in order to determine the signal decay: where x is the distance below the water surface and bBS describes the backscatter. The function approximation delivers an estimate of the volume exponential decay coefficient k.

SVB algorithm
The SVB algorithm extracts the echoes of water surface, water column and water bottom (Schwarz et al., 2019). It is based on the exponential decomposition proposed in Schwarz et al. (2017), which is a waveform decomposition method using a physically motivated model of the dBCS and a record of the system waveform. The dBCS model is composed of two boxcar shaped function segments for water surface and bottom as well as one exponential function segment for the water column ( where tn are the sampling time instances of the sample set of size N . The SVB algorithm delivers the parameters ρ of the dBCS model. The exponential decay coefficient of the function segment for the water column is used as a measure of water turbidity. The approach has the advantage that neither the backscatter at the water surface nor that at the water bottom distort the determination of the decay coefficient. Furthermore, the influence of the system waveform on the turbidity determination is eliminated. The International Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences, Volume XLIII-B2-2021 XXIV ISPRS Congress (2021 edition) time intensity Figure 3. Schematic representation of SVB algorithm; green: dBCS model; red: model of the received signal; black: received signal.

DATA
For practical evaluation we used a LiDAR bathymetry dataset of the Elbe River, provided by German Federal Waterways and shipping Office Dresden (WSA) and the German Federal Institute of Hydrology (BfG). The Elbe River is a shallow inland water regulated by numerous groynes in the study area. A 156 m long and 154 m wide river section was selected for the investigations, with a water depth of about 2.1 m in the shipping channel. In the selected test area, there is also a groyne which influences the flow velocity. The ALB data was acquired in spring 2015 with a RIEGL VQ-880-G (Riegl, 2021) using a palmer scan pattern with 20 • incidence angle in a flying height of 380 m above ground. The point density is about 55 points/m 2 caused by overlapping flight strips.
The two water turbidity estimation methods are applied both to the original measurement data (individual full-wavforms) and to further processed data (stacked full-waveforms and ortho full-wavforms). The data processing and characteristics of the full-wavform types are described in 3.1 to 3.3.

Individual full-waveforms
The individual full-waveforms are the unprocessed recorded full-waveforms consiting of 60 to 200 samples (intensity values). The constant sampling interval is 0.575 ns, corresponding to a sampling distance of 6.45 cm under water. The fullwaveform signals result from the propagation of the laser pulse in the water column, which is oblique. Therefore, turbidity parameters derived from those signals do not correspond to turbidity in the vertical direction. Figure 4 shows a typical individual full-waveform, which is affected by some noise.

Stacked full-waveforms
Stacked full-waveforms are obtained by aggregating closely adjacent individual full-waveforms with similar properties with regard to the water depth. As presented in Mader et al. (2019Mader et al. ( , 2021, data processing starts with the partitioning of the measured data into a regular grid. The grid cell size depends on the data density and the ground topography characteristics (strong or quick changes in terrain slope). A grid cell size of 2 m × 2 m was used for the present data set, resulting in 109 to 144 waveforms per cell.
Combining the individual full-waveforms in a grid cell requires the correct alignment to each other with respect to the height coordinates. Since information on the position of pulse emission, laser beam direction and the laser pulse travel time was not available for our data, the vertical pulse alignment was realized using the maxima of the water surface echoes in the fullwaveforms. This approach is based on the assumption that the water surface approximately corresponds to a horizontal plane at the time of measurement. Subsequently, all measured fullwaveform signals of the respective grid cell are summed up resulting in a stacked full-waveform ( Figure 5).  Figure 6 shows a typical stacked full-waveform, which has a significantly lower signal noise than an individual fullwaveform. Due to the principle of the full-waveform stacking approach, full-wavforms of the forward look and the backward look of the palmer scan are jointly evaluated. Therefore, the stacked full-waveform signals result from the averaging of different oblique laser beam directions. The derivation of turbidity parameters from stacked full-waveforms is thus based on the assumption that the turbidity within the grid cell is independent of the angle of incidence of the laser pulse.

Ortho full-waveforms
An alternative approach of the combined evaluation of fullwaveforms is the generation of ortho full-waveforms (Pan et al., 2016). For this purpose, the individual full waveforms are integrated into a voxel space. The basic idea is the voxelization of the irregular point cloud, which consists of the georeferenced samples of all individual full-waveforms, i.e. the approach is based on a narrow laser beam model . The voxel attributes result from the mean intensity of the samples included in the voxel. In order to generate ortho full-waveforms, we analyze the vertically superimposed voxels representing the vertical water column. More precisely, an artificial full-waveform is assembled from the voxel attributes of each voxel space column. Figure 7 visualizes the principle of ortho full-waveforms generation. For the present data set the voxel size was set non-cubic to 2 m × 2 m × 0.1 m. The vertical extent of the voxels was chosen to be larger than the sampling interval of the individual fullwaveforms to ensure that a waveform is represented at least once in each voxel along the voxel space column. The number of individual full-waveforms used to generate one ortho fullwaveforms is not clearly defined. Due to the oblique beam path, an individual full-waveform can intersect several voxel columns. Between 174 to 446 individual waveforms contribute to a ortho full-waveform. The ortho full-waveforms initially have a sampling distance of 0.1 m, corresponding to the vertical extension of the individual voxels. The application of the SVB algorithm requires re-sampling in order to adapt the sampling interval to that of the system waveform. Figure 8 shows a typical ortho full-waveform, which is characterized by an improved signal-to-noise ratio compared to individual full-waveforms. The use of ortho full-waveforms enables the strictly vertical derivation of turbidity information.

RESULTS AND DISCUSSION
The two methods for deriving turbidity parameters were applied to all three waveform types. The results are presented and discussed in Section 4.1. Section 4.2 presents the results of the quantitative evaluation. Figure 9 shows the derivation of the turbidity parameters with exponential function approximation (top) and SVB algorithm (middle and bottom). The high signal noise of the individual full-waveform influences both the quality of the exponential function approximation and the determination of the volume backscatter with the SVB algorithm. Using stacked full-waveforms and ortho full-waveforms comes with the advantage of lower signal noise, providing a more reliable determination of turbidity.

Turbidity estimation
Strictly speaking, the SVB algorithm based on the exponential decomposition is only suitable for the application to individual full-waveforms that have a direct physical relation between received waveform and system waveform. The application of the algorithm to further processed waveform data is based on the assumption that the influence of the system waveform on the received waveform is not changed by stacking or voxelization.
The exponential function approximation works correctly for the majority of the individual, stacked, and ortho full-waveforms. However, there may be problems when the water column has not been reliably extracted. This concerns especially the water bottom echo. For the ortho full-waveforms, there is currently no control and correction of a possibly incorrectly determined water bottom echo, as it is done for the individual full-waveforms and stacked full-waveforms.
With the SVB algorithm, problems occur more frequently. In many cases, the echoes from the water surface, water column and water bottom are not extracted correctly. As a result, the estimated model of the received signal does not approximate the received signal well. This affects the derivation of the exponential decay coefficient of the volume backscatter. We expect that the overall performance of the SVB algorithm can be improved by optimizing the determination of approximate values. Further investigations are necessary to adapt the algorithm to the available data.
In addition, problems arise in the shore area where the waveforms deviate from the characteristic signal shape due to the shallow water depth. In very shallow water, the echoes from the water surface, water column, and bottom superimpose strongly (Wagner et al., 2006;Ullrich and Pfennigbauer, 2011). In some cases, the echo from the water bottom even has a larger amplitude than the echo from the water surface. In these special cases, the derivation of the decay coefficient fails. However, the developed method is able to detect such cases automatically. Figure 10 shows the results of the turbidity estimation in the study area in a color-coded visualization. The XY-coordinates correspond to the water surface points, the decay coefficient k is color-coded from blue (low turbidity → clear water) to brown (high turbidity → turbid water). For the exponential function approximation, a small value for k indicates low turbidity and a large value indicates high turbidity. In the SVB method it is the opposite. The stronger the exponential decay in the received model, the weaker the decay in the dBCS (see figure 9). A small value for k indicates a high turbidity and vice versa. For this reason, the color bar is inverted in the SVB algorithm.
Due to the correlation of the decay parameter with other water column parameters of the dBCS model, further investigations are necessary.
The evaluation of individual full-waveforms achieves a significantly higher spatial resolution than the evaluation of stacked full-waveforms and ortho full-waveforms. Compared to common in-situ measurement methods, which only allow a pointwise detection of turbidity parameters, the resolution of the further processed full-waveform types is still very high. The results show some local turbidity variations within the river. Noticeable are the clearer areas in the middle as well as near the groyne. Basically, the determined turbidity values are similar for the different full-waveform types as well as for the two turbidity estimation methods.
With the exponential function approximation applied to individual full-waveforms, the values for k range from 1.3 m −1 to 3.4 m −1 . As expected, the results for stacked full-waveforms and ortho full-waveforms show a lower spread. They have a similar range from 1.8 m −1 to 3.2 m −1 . The small differences can be explained by the vertical derivation of the turbidity parameters for the ortho full-waveforms.
The SVB algorithm generally provides less plausible results, characterized by a strong noise. For the stacked full-waveforms and the ortho full-waveforms, the turbidity values are in the same range as for the exponential function approximation. The values for the individual full waveforms vary between 0.5 m −1 and 3.5 m −1 . The local turbidity variations are basically comparable to the results obtained with the exponential function approximation.

Evaluation
The area-wide collection of turbidity ground truth data requires a dense grid of in-situ measurements. The acquisition of these data is very time-consuming and cost-intensive. Furthermore, it must take place at the same time as the flight. In addition, conventional turbidity measurement methods determine different optical properties of the water whose physical relationship to the turbidity parameters obtained from ALB data is not trivial. In our study no turbidity ground truth data are available. Therefore, the comparison of the turbidity estimation methods applied on the different waveform types is based on a quantitative evaluation of the derived turbidity parameters.
For evaluation purposes, the correlations between adjacent turbidity parameters are investigated. Since the turbidity in the medium water will not change abruptly, neighboring fullwaveforms should deliver a similar turbidity. To investigate the correlation we form pairs of values consisting of turbidity parameter and neighboring turbidity parameter. For individual full-waveforms we use the nearest neighbor, while for stacked full-waveforms and ortho full-waveforms all 8 neighbors in the grid are considered separately. Each pair of values is taken into account only once. Figure 11 and Table 1 show the results of the correlation analysis. Please note that a complete agreement between the turbidity parameters derived from neighboring full-waveforms cannot be achieved since a small change in turbidity parameters is expected due to local turbidity variations. Looking at Figure 11 one can notice that the majority of the data points in the scatter plots generally form a cluster. In addition, there are a number of data points that deviate more or less strongly from the 1:1 line. In principle, the cluster is more compact with the exponential function approximation and the spread of the data points is smaller.
The correlation coefficient values (Table 1) range between 0.12 and 0.58. The correlations for the SVB algorithm are relatively poor for individual and ortho full-waveforms. The SVB algorithm applied on stacked full-waveforms results in similar correlation like the exponential function approximation. The highest correlation occurs with exponential function approximation applied on stacked full-waveforms and ortho fullwaveforms. This indicates that the neighboring turbidity values are very similar proving the plausibility of the results.

CONCLUSION
Airborne LiDAR bathymetry represents an interesting option for the derivation of area-wide turbidity information for applications in both limnological studies of inland waters as well as oceanographic studies. Analyzing the voulme backscatter of the ALB full-waveforms allows for the determination of the signal decay, which represents a measure for water turbidity. Compared to conventional point-wise turbidity estimation methods, a higher degree of automation and an area-wise data capture at a much higher spatial resolution can be achieved.
The paper compares two different methods for the estimation of the exponential decay coefficient k, namely the exponential function approximation and the SVB algorithm. Both approaches are applied to three different full-waveform types (individual, stacked, and ortho full-waveform). Due to the lack of ground truth data, an external validation of the methods was not possible yet. However, a quantitative evaluation of the results showed that the exponential function approximation applied on stacked full-waveforms and ortho full-waveforms and the SVB algorith apllied on stacked full-waveforms performs best for the derivation of turbidity parameters. Since the use of ortho fullwaveforms enables the vertical derivation of turbidity parameters we favor the exponential function approximation applied on ortho full-waveforms.
The aim of future work is to confirm the validity of the turbidity determination in further experiments by acquiring terrestrial reference data. In addition, the extraction of the water column from the ortho full-waveforms will be improved by developing volumetric approaches to derive water bottom information. Moreover, test areas with a water depth larger than 2 m will be investigated. There, the volume backscatter consists of significantly more samples, enabling a more reliable derivation of the turbidity parameters. At the moment the exponential function approximation is limited to the derivation of one turbidity parameter per waveform assuming a uniform turbidity distribution within the vertical water column. Future work will focus on vertical variations in turbidity by examining the volume backscatter segment-wise.