DETERMINATION OF 3D WATER TURBIDITY PARAMETER FIELDS FROM LIDAR BATHYMETRY DATA BY VOLUMETRIC DATA ANALYSIS

: Accurate information on turbidity in water bodies is relevant to numerous limnological and oceanological issues. However, the collection of turbidity parameters using conventional in-situ measurement methods is time-consuming and cost-intensive and therefore usually limited to very small study areas. The use of airborne LiDAR bathymetry data is a promising alternative. However, existing methods for deriving turbidity parameters from airborne LiDAR bathymetry data are limited to the determination of one single turbidity parameter per water column element. The paper presents a novel approach that overcomes the existing limitations enables the determination of 3D water turbidity ﬁelds. By volumetric data analysis, the vertical turbidity stratiﬁcation in the water body can be determined. For validation purposes, the approach was applied to synthetic measurement data generated in a simulation as well as a real measurement data set of a shallow coastal water.


INTRODUCTION
Airborne LiDAR bathymetry (ALB) is an efficient technique for surveying shallow water bodies. Water turbidity (due to sediment load, algae, pollutants, etc.) represents a crucial influence parameter for the depth performance of the measurement method. Turbidity causes a continuous decrease of the signal intensity. Conversely, an analysis of the signal decay can therefore provide quantitative information about the turbidity of the water column penetrated by the laser pulse.
Turbidity is not constant for the whole water body but can show local variations in horizontal and vertical direction. The degree of turbidity depends on the properties and distribution of the scattering matter, which is closely related to the ecological processes and sediment transportation in the water body. Scattering particles moving in the water, for example, are an indicator of the dynamics of the upper water layer (Churnside, 2014). Furthermore, nutrient input, freshwater supply and temperature changes can create thin layers with a high concentration of organisms or particles, which result in higher water turbidity (Sullivan et al., 2012). Therefore, information on the vertical turbidity stratification allows conclusions on ecological processes and the properties of the water body.
In Richter et al. (2021b) we present an approach to derive turbidity parameters from LiDAR bathymetry data that builds on previous work published in Richter et al. (2017). The method is based on the analysis of the water column part of the received waveform signals by approximating an exponential function. Herein the exponential decay coefficient represents a quantitative integral non-metric measure of water turbidity. Compared * Corresponding author to conventional point-wise measurement methods for determining water turbidity (e.g. measuring of visibility depth with a Secchi-disk, measuring of down welling irradiance with submersible quantum sensors, measuring of a beam attenuation coefficient with transmissometers), the approach has the advantage of providing full field area data. However, only one turbidity value is determined for each full-waveform signal, i.e. information on the vertical turbidity stratification is not available. This paper will present a first approach for the derivation of depth-resolved turbidity parameters. The aim is to move from an area-based determination of water turbidity to a volumetric representation in order to additionally record the depth stratification of the water turbidity. The point density of modern 'small footprint' LiDAR bathymetry systems is used to determine spatially resolved decay parameter fields in a voxel space representation of the LiDAR bathymetry data. We expect that this approach can achieve significant results at least at larger penetration depths and pronounced stratification of the water turbidity. The methodology of the approach is described in Section 2. Section 3 presents the data sets available for testing the method. The results of the practical investigations are summarized in Section 4. Section 5 deals with the problems of validation. The paper finishes with a short conclusion in Section 6.

METHODS
The basis for the determination of depth-resolved turbidity parameter fields is the generation of a local voxel space representation from neighboring digitized laser pulse echoes. For this purpose, the individual full-waveforms are integrated into a voxel space by voxelising the point cloud of the georeferenced samples of all individual full-waveforms. The approach is therefore based on a narrow laser beam model (Richter et al., 2021a). The voxel attributes result from the mean intensity of the samples included in the voxel. Subsequently, ortho-full waveforms are derived from the voxel space representation by analyzing the vertically superimposed voxels representing the vertical water column (Pan et al., 2016;Richter et al., 2021b;Mader et al., 2022). More precisely, a vertical pseudo fullwaveform is assembled from the voxel attributes of each voxel space column. Figure 1 visualizes the principle of ortho fullwaveforms generation.
In the next step, the water column part (volume backscatter) of the ortho full-waveform is extracted. It 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. To reliably eliminate the influence of the water surface echo and water bottom echo, the water column is shortened by 10 samples at both ends.
Finally, depth-resolved turbidity parameters are derived for each extracted volume backscatter by fitting multiple exponential segments into the volume backscatter. The principle is shown schematically in Figure 2. The number of required segments is derived from the shape of the volume backscatter in an iterative process. Initially, starting from the boundary of the water surface echo, only the fist part of the volume backscatter is considered. Theoretically, 3 samples would be sufficient for the function approximation. However, to obtain a more reliable result that is less affected by noise, 5 samples are chosen. An exponential function is fitted into this part: 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. The goodness of the fit can be assessed via the Root Mean Square (RMS).
Subsequently, further samples are added step by step and the function approximation is repeated until the end of the volume backscatter is reached. Then the RMS values of all function approximations are analyzed. The following cases may occur:  1. The turbidity is constant in the considered part of the water column. The data is well described by the functional model and the RMS is small.
2. The considered part of the water column contains few samples of a turbidity variation. The data is worse described by the functional model and the RMS deteriorates.
3. The considered part of the water column contains many samples of a turbidity variation. The data is poorly described by the functional model and the RMS is large. Figure 3 shows the situation schematically. An empirically determined threshold value can be used to specify the number of considered samples above which the RMS value increases. A new exponential segment is added at this location. The whole process is repeated until all exponential segments have been found.

DATA
For evaluation purposes, synthetic ortho full-waveforms generated in a simulation (Section 3.1) as well as real measurement data from a shallow coastal water (Section 3.2) are available. This section describes the data basis in detail.

Simulation
The simulation generates ortho full-waveforms of an 8 m deep water body with three layers of different turbidity. The position, extent and turbidity of the layers are chosen randomly within certain limits. The position of the second turbidity layer is selected so that the upper and lower limits are at least 1 m from the water surface or the bottom of the water body. The thickness of the second turbidity layer varies between 1 m and 3 m. In order to generate the synthetic ortho full-waveforms, a model of the differential backscatter cross section (dBCS) according to Schwarz et al. (2017) is set up first. The dBCS is composed of several exponential segments of the form: For γ = 0 a boxcar shaped segment results, which can be used for the modeling of water surface and water bottom.
In the present scenario, the dBCS is composed of the following segments: one boxcar shaped function for the water surface, three exponentials for the different layers of the water column, and one boxcar shaped function for the water bottom. The positions τi and widths Ti of the segments result from the selected water depth and the position of the second turbidity layer. The width of the water surface segment and the water bottom segment is set to T = 4. Figure 4 shows a resulting dBCS.
The ortho full-waveform results from the convolution of the modeled dBCS and a real system waveform from a Riegl VQ-880 ALB system ( Figure 5). The convoluted signal is resampled at an interval of 0.575 ns which corresponds to a distance of 6.47 cm in water at a speed of light of 225000000 m s . To represent the character of real measurement data, we add white Gaussian noise to the discrete samples of the ortho full-waveform as a simple noise model for modeling the noise behaviour of a receiver. The amplitude of the noise can be modeled with a normal probability distribution. Figure 6 shows the resulting ortho full-waveform. A total of 1000 simulations were performed in this way.

Real measurement data
The measurement data were acquired in June 2020 with the ALB system Leica Chiroptera 4X in a shallow coastal water area. The study area is located east of the island of Sylt in the North Sea. The water depth in the study area is 3 m to 5 m. Figure 7 shows an ortho full-waveform generated from the volumetric representation of the raw data. In the study area, 300 ortho full-waveforms were created for further investigation.   . Ortho full-waveform (black) resulting from the volumetric evaluation of real measurement data (Richter et al., 2021b;Mader et al., 2022) and water column (blue) extracted with full-waveform processing approach presented in Mader et al. (2019Mader et al. ( , 2021.  Figure 8 shows the result of the segmentwise exponential function approximation applied on the simulated ortho fullwaveform from Figure 6. Three exponential segments were fitted into the water column part of the signal. The positions of the two turbidity changes result from the transition points between the exponential segments. The distances between water surface and transition points can be converted to a metric value using the known sample distance of the ortho full-waveform. The same applies to the thickness of the turbidity layers, which results from the widths of the exponential segments. The determined parameters are summarized in Table 1.

RESULTS
To evaluate the result, Figure 8 also shows the actual position and thickness of the turbidity layers (blue rectangles) as known from the simulation. The comparison to the second exponential segment reveals an offset of 5.175 ns which corresponds to a distance of 33.5 cm. A significant change in the RMS value can only be detected when a sufficiently large number of samples of the next turbidity layer is included in the function approximation. Therefore, several samples are always assigned to the previous segment, causing a certain backlag which will be addressed in future work.
The result of the segment-wise exponential function approximation applied to the real data ortho full-waveform from Figure  7 is shown in Figure 9. In contrast to the simulated data, only 2 exponential segments were fitted here because the water body in the study area has only one turbidity variation in the depth direction. The position of the turbidity change, the thickness of the turbidity layers, and the decay parameters are summarized in Table 2.

VALIDATION
For the simulated ortho full-waveforms, location, thickness, and decay parameter of the turbidity layers are precisely known 1.00 decay parameter k1 0.042 decay parameter k2 0.064 Table 2. Parameters estimated from real ortho full-waveform.
from the simulation and can be used for validation. For this purpose, the deviations between reference values and actual values are statistically evaluated.
The histograms of the location deviations in Figure 10 (a) and (b) show that the position of the turbidity change is systematically determined about 10 cm to 20 cm too deep. This is due to the fact that some samples of the turbidity layer were assigned to the previous layer. 50.1 % of the values deviate less than ±20 cm from the reference location.
The thickness of the turbidity layers tends to be determined too low as shown in the histograms in Figure 10 (c)-(e). In the first and third segment, this is due to the extraction of the water column at a fixed distance from the detected water surface and water bottom echo. In order to remove both echoes completely, this distance is chosen quite generously. This leads to the fact that the extracted water column is generally too short. Nevertheless, the second segment is not affected by this fact. As shown in the histogram in Figure 10 (d), however, the thickness is determined about 10 cm to 20 cm too small. 52.0 % of the values deviate less than ±20 cm from the reference thickness.
The deviation of the decay parameter is given in percentage of the reference value. The decay parameters in the first and second segment show small deviations of ±5 %, which have a normal distribution (Figure 10 (f) and (g)). In the third exponential segment, the decay parameter tends to be determined too small. 73.9 % of the values, however, deviate less than ±5 % from the reference value.
With real measurement data, validation is much more difficult. The collection of depth resolved turbidity ground truth data requires a large number of in-situ measurements, which is very time-consuming and cost-intensive. In our study no turbidity ground truth data were available. To evaluate the methodology, we therefore examine the consistency of the results in a close neighborhood. Since the turbidity in the medium water will not change abruptly, neighboring ortho full-waveforms should deliver similar positions of turbidity change and decay parameters. No lateral turbidity variations are expected in the smallscale study area. maximum at 2.2 m. Since the beginning of the water column is at a constant distance from the detected water surface echo, the histogram of the thickness of layer 1 in Figure 11 (b) shows an identical characteristic. The maximum is located at 2.2 m. The thickness of layer 2 varies between 2 m and 3 m ( Figure  11 (c)). Obviously, the position of the turbidity change and the thickness of the turbidity layers cannot be determined with great certainty. The reason for this is the weak expression of the turbidity variation. The histograms of the decay parameters in Figure 11 (b) and (c) show that there is little difference between the turbidity of the two layers.

CONCLUSION
The paper presents a method for the determination of 3D water turbidity parameter fields frim LiDAR bathymetry data. The basic idea is the generation of vertical ortho full-waveforms, which are derived from a volumetric representation of the measurement data. Depth-resolved water turbidity is determined by analyzing the signal decay in the ortho full-waveforms. For this purpose, multiple exponential segments are fitted into the water column part of the signal. The number of required segments and the boundaries between the segments are derived from the ortho full-waveform signal itself.
Validation of the developed methodology with simulated data provides promising results. When tested on real measurement data, the method reached its limits. The turbidity stratification present in the study area was very weak and therefore could not be reliably detected. However, it is expected that the method can be successfully applied to real measured data if the turbidity variations are pronounced enough.
In future work, therefore, measurement data should be acquired from a water body with more prominent turbidity stratification. The simultaneous acquisition of ground truth data with conventional in-situ turbidity measurement methods will allow a comprehensive validation of the results.