NUMERICAL SIMULATION AND EXPERIMENTAL VALIDATION OF WAVE PATTERN INDUCED COORDINATE ERRORS IN AIRBORNE LIDAR BATHYMETRY

Airborne LiDAR bathymetry (ALB) requires a refraction correction on the basis of Snell’s law at the air-water interface and a speedof-light correction to be applied on the raw laser data in order to achieve a geometric accurate representation of the water bottom. Strictly speaking, this requires exact knowledge about the local water surface inclination. If this information is not available, certain simplifications have to be introduced in correction methods. Common correction methods assume either a horizontal or a locally tilted planar water surface as well as an infinitesimally small thin laser ray, thus neglecting effects caused by the finite laser pulse diameter penetrating a curved surface. In our simulation approach, the refraction of finite diameter laser pulses passing the air/water interface is modeled differentially in a strict manner. The simulation tool is able to predict wave induced coordinate errors which have to be expected due to the neglections made in common refraction correction methods. Moreover, wave pattern dependent correction terms were be derived from systematic portions of the errors revealed by the simulations. The goal of this paper is to experimentally validate the coordinate errors predicted by the simulation tool. For that purpose, airborne laser bathymetry data of a 12 by 50 meter open air wave pool were processed, and the results were compared to reference data of the empty pool acquired by terrestrial laser scanning. The comparison showed that the effects predicted in the numerical simulation are confirmed by the experimental validation.


INTRODUCTION
Geometric modeling in airborne LiDAR bathymetry is significantly more complex than in conventional laser scanning since the laser pulse is passing two different media.Due to the different refraction indices the direction and velocity of the laser pulse propagation is changing at the interface surface between the two media.Additionally, the laser ray is subjected to influences such as beam spreading due to dispersion at small sedimentary particles and organic materials and the diffuse reflection at the water bottom, both contribute to a different light path towards the sensor.Figure 1 shows the refraction of the incident laser pulse based on the local wave-induced water surface.In order to achieve a geometric accurate representation of the water bottom a run-time correction and a refraction correction on the basis of Snell's law has to be applied on the raw laser data.The refraction correction requires exact information on the local water surface inclination.As this information is usually not available, existing correction methods introduce certain simplifications.
The most common correction method assumes a horizontal and planar water surface at which the laser beam is refracted (Fig. 1, purple).Even small deviations from the planarity, already caused by moderate sea swell, can lead to a significant lateral displacement dXY hz and height displacement dZ hz .More complex correction methods try to consider the actual water surface geometry.For this purpose, some system providers apply a local water surface tilt (Fig. 1, red), e.g. based on the intersection of the incident laser ray with a triangular mesh of detected water surface points * Corresponding author in their bathymetry data processing methods (Ullrich and Pfennigbauer, 2011).If such a linear local tilt has been considered in data processing, the effects of wave patterns on coordinate accuracy will be reduced and the coordinate displacements dXY tilt and dZ tilt will be smaller.
Beside the simplifications concerning the water surface geometry, both correction approaches are limited by the fact that the laser ray is considered to be an infinitesimal small line only.Effects caused by a finite diameter laser pulse penetrating a curved surface are neglected.In contrast, our simulation approach investigates the effect of waves patterns on LiDAR bathymetry water body bottom coordinates under strict consideration of refraction effects (Westfeld et al., 2017).For this purpose, the refraction of finite diameter laser pulses passing the air/water interface is modeled differentially in a strict manner.The developed models can be used to correct each individual laser pulse.Therefore, the reconstruction of the actual water surface, for instance from dense laser scanner points, is required.If this information is not available, the simulation results can also be used to derive wave pattern dependent correction terms which cover systematic errors.
In our previous simulation studies we examined the effect of wave patterns on refraction and subsequently on coordinate accuracy for typical ocean wave patterns (Westfeld et al., 2017) as well as riverine wave patterns (Westfeld et al., 2016).It has been shown that, depending on sea swell and laser footprint size, the effect on lateral bottom point displacement dXY hz can amount to several decimeters, in some cases even meters in the planimetry coordinates of underwater points.Furthermore, height displacements dZ hz in decimeter range have to be taken into account.This contribution complements the numerical simulation by an experimental validation of wave pattern induced coordinate errors for a real world scenario.The aim of the paper is to examine how the coordinate errors predicted in our simulation correspond to the errors derived from real measurement data acquired from a 12 by 50 meter open air wave pool.For this purpose, we apply the most simple refraction correction method assuming a horizontal water surface as well as the more complex refraction correction method with local surface tilt on the raw measurement data.We compare the refraction corrected data with terrestrial reference data to assess the coordinate errors remaining in the data after conventional refraction correction.The depth coordinates displacement is derived from the data on the pool bottom whereas planimetric coordinate displacements can be determined from points on the pool wall.In addition to the analysis of the measurement data we use our simulation to predict systematic coordinate errors with respect to both refraction correction methods mentioned above.For this purpose, we expand our modeling approach to the local water surface tilt based on our previous findings referring to a horizontal water surface.By choosing suitable simulation parameters we can reproduce the wave pattern like it is actually present in the measurement data.
The paper is structured in the following way: Section 2 briefly describes the experimental setup and the acquired reference and measurement data.The methods for numerical simulation and experimental validation are given in section 3 and 4. In section 5, the results are presented and discussed whereas section 6 summarizes the work and addresses future tasks.

EXPERIMENTAL SETUP
For the experimental investigations we chose an open air wave pool which has a total length of 50 m and a width of 12 m.With its wave generation possibilities as well as with its regular geometry and good reflection properties, the pool seems well suited for a validation study.The bottom of the pool slopes from 0.3 m to a maximum water depth of 1.5 m (at horizontal water surface).The wave machine produces regular waves with amplitudes of about 0.5 m, which are characterized by braking crests and whitecaps especially in shallow water.The experimental setup offers controlled examination conditions by reliably producing waves with known parameters, low water turbidity and a fixed precisely measurable water bottom geometry.
The reference data was collected by terrestrial laser scanning while the wave pool was empty.To achieve a dense representation of the water bottom geometry the pool was scanned with a Riegl LMS-Z420i from five different positions.Figure 2 shows the data acquisition as well as the resulting reference point cloud.
The airborne survey campaign was carried out in the filled state under wavy as well as smooth water surface conditions (fig.3).The LiDAR bathymetry data was acquired with a RIEGL VQ 820G LiDAR system in different flying heights (500 m, 600 m, 700 m) and flight directions (in direction of wave propagation and across).The laser beam divergence of 1 mrad results in a laser footprint with a diameter of approximately 0.5 m to 0.7 m at the water surface.Figure 3 shows a profile of the ALB point cloud under wavy conditions.

NUMERICAL SIMULATION
The numerical simulation of wave pattern induced effects on refraction and thus on the planimetric and depth coordinates of water bottom points requires water surface modeling, bottom surface modeling and ray path modeling.
For the water surface modeling we used Tessendorf's oceanographic statistics based surface wave model for ocean waves (Tessendorf, 2001).The model treats each wave height as a random variable of its planimetric position at a given time t.Based on Tessendorf's model a height field is generated by means of Fast Fourier transform (FFT).The height field represents a realistic ocean surface in the form of a dense regular grid.The charac- teristics of the height field can be influenced by parameters of the Fourier grid (width, height, mesh size) as well as wind conditions (wind speed, wind direction).In order to achieve comparability to the experimental validation, we aim at reproducing the wave pattern like it is actually present in the measurement data.For this purpose, we analyzed the wave pattern represented by the measured water surface points to derive its amplitude and wavelength.The wave amplitude refers to the vertical distance from mean level to crest and the wave length specifies the horizontal distance from crest to crest.Subsequently, we choose suitable simulation parameters to obtain a water surface model with similar properties.
The bottom surface modeling is focused on the plane characteristic of the actual pool bottom.We deliberately omit the slope down, since the predicted measurement errors will be specified in percent of the water depth.Therefore, the simulated water bottom is generated as a horizontal plane surface.

EXPERIMENTAL VALIDATION
The experimental validation is based on the LiDAR bathymetry data as well as the terrestrial laser scanner data, which serve as reference for the following tests.The bathymetry data, acquired in the airborne survey campaign, is provided as uncorrected 3D point cloud, i. e. no refraction correction and run time correction was applied on the raw data set.An accurate time stamp is available for each 3D point in addition to the classification in water surface and water bottom points.Furthermore, information on the sensors trajectory and manufacturer specifications regarding the refractive indices of air and water are accessible.
Based on these information we apply a runtime correction and the simple refraction correction method assuming a horizontal water surface as well as the more complex refraction correction method The International Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences, Volume XLII-2, 2018 ISPRS TC II Mid-term Symposium "Towards Photogrammetry 2020", 4-7 June 2018, Riva del Garda, Italy with local surface tilt on the raw measurement data.The assumption of a horizontal water surface can be achieved by either a constant mean water level height over the entire area of investigation or by local horizontally oriented water surface elements at different heights provided by the water surface points.The following cases are analyzed: 1. horizontally oriented water surface elements and constant mean water level height (M1) 2. horizontally oriented water surface elements and different local heights (M2) 3. tilted water surface elements and different local heights (M3) In the first case, the water surface is represented by a plane, whose height is extracted from the measurement data acquired under smooth water surface conditions.In the two other cases we generate a mesh including the original water surface points with linear interpolation methods.The water surface is represented by a triangulation of the mesh points.To realize the refraction correction we estimate the direction vector between each water bottom point and the corresponding trajectory point.Subsequently, the direction vector is intersected with the water surface.Please note, that the laser pulse is treated as an infinitesimal small line here.In case of the meshed water surface the local height is derived by linear interpolation between the vertices of the intersected triangle.Finally, the direction of the laser ray is corrected applying Snell's law.Considering the reduced velocity of light in water results in the corrected water bottom point coordinates.To evaluate the correction results we analyze the differences between the coordinates obtained from the different correction methods.Furthermore, we compare the refraction corrected data with the terrestrial reference data to assess the coordinate errors remaining in the data after conventional refraction correction.The depth coordinates displacement is derived from the data on the pool bottom whereas planimetric coordinate displacements can be determined from points on the pool wall.In order to reference ALS and TLS point clouds we perform a rough registration with three homologous points in both point clouds followed by an ICPbased fine registration.The registration accuracy is in the range of several centimeters.

Numerical Simulation
In order to simulate the wave pattern, like it is actual present in the measurement data, we analyze the measured water surface points.For this purpose, we aggregate 50 cm wide sections of the water surface point cloud to profiles.Afterwards, we derive wave parameters by fitting a spline function into each profile.The maximum amplitude and the wave length arise from the local minima and maxima, representing wave crests and troughs.Figure 4 (a) shows a typical profile with a maximal wave amplitude of absolutely 0.96 m (max.wave crest height 0.50 m, min.wave trough height −0.46 m).The maximal wave length is 8.00 m.The origin of the waves is on the right side, whereby the water depth grows with increasing X-coordinates.The mean water level is visualized as a horizontal line.Using appropriate simulation parameters, the actual wave pattern is reproduced as close as possible.Figure 4  wave troughs is −0.31 m resulting in an amplitude of absolutely 0.77 m.The maximal wave length is 10 m.
The simulation is run for 1000 consecutive epochs, taking into consideration the refraction at the local wave-induced water surface as well as the refraction at the horizontal or locally tilted water surface.The point density of the locally tilted water surface defines the representation accuracy of the triangulation.Due to the inhomogeneous distribution of the water surface points in the measurement data we consider two cases with 1 point per square meter (p/m 2 ) as well as 10 p/m 2 .
The resulting coordinate displacements are presented in figure 5 and table 1.As the effect of wave patterns on refraction increases linearly with water depth, all results are presented in percentage of the water depth.The coordinate displacements consist of both a lateral component dXY (red curve in fig.5, row 1-3 in table 1) and a depth component dZ (blue curve in fig.5, row 4-6 in table 1).
The root mean square error (RMSE) of the lateral coordinate displacement at a flying height of 500 m is 1.40 % (max.3.12 %) of the water depth for the horizontal water surface.The locally tilted water surface with a point density of 1 p/m 2 results in a RMSE of 1.02 % (max.2.45 %).Assuming a water depth of 1.6 m, a RMSE of 1.6 cm (max.3.9 cm) has to be expected in areas with low point density.The RMSE is reduced to 0.23 % (max.0.64 %) corresponding to 0.4 cm (max.1.0 cm) if the tilted water surface is represented by 10 p/m 2 .In summary, the lateral coordinate displacements decrease with increasing complexity of the water surface.The same applies to the other two flying heights, whereby the coordinate errors in planimetry decrease with increasing laser beam footprint.
In general, the depth component of the coordinate displacement is The less prominent than the lateral component.For a flying height of 500 m, the RMS depth errors vary in a range of 0.12 % to 0.28 % (0.2 cm to 0.4 cm) depending on the chosen correction method.
The variations between the different flying heights are less distinctive.Furthermore, the water surface complexity does not necessarily affect the depth errors (cf.hz600/t1600).

Experimental Validation
In contrast to the simulation approach, the local wave-induced water surface (Fig. 1, blue) is unknown in the experimental validation, inhibiting the calculation of dXYhz, dZhz, dXYtilt and dZtilt.Therefore, the discrepancies between the coordinates obtained from the different correction methods are analyzed.The results presented in Table 2 show that the planimetric coordinates are mainly influenced by the local surface tilt.At a flying height of 500 m the consideration of the local height (M1 -M2) reveals a coordinate difference of 2.29 %.Taking into account the local tilt of the water surface results in a difference of 6.67 % (M1 -M3) respectively 6.17 % (M2 -M3  2. Discrepancies between the coordinates obtained from the different correction methods Mi (M1 -M2, M1 -M3, M2 -M3) in percent of the water depth for different flying heights (500 m, 600 m, 700 m).
coordinate differences of 1.30 % for the comparison of method M2 and M3 show that the local surface tilt is less important for the depth coordinate whereas the local height is more essential (M1 -M2 and M1 -M3).The results are similar for different flying heights.
The comparison of refraction corrected data and terrestrial reference data is based on the investigation of deviations between ALS and TLS point cloud at the pool bottom (depth displacement) and at the pool wall (planimetric displacement).Due to the oblique incidence angle in combination with the highly reflecting material, the point coordinates at the pool wall are affected by noise which superimposes the geometric effects.Therefore, only a small number of points at the concrete base of the water slide are available for the planimetry effect evaluation.The relevant region is marked in figure 6.The deviations initially contain systematic errors due to the limited registration accuracy as well as effects of the beam divergence and incidence angle.In order to achieve the random part of the deviations, the systematic errors have to be eliminated.For that purpose, we estimate the systematic errors based on the results of the third refraction correction method M3, which provides the best available correction results (cf.table 2).point clouds and TLS reference point cloud is limited by the small number of usable ALS points solely available for the 500 m flight strips.Figure 9 presents the results for the different refraction correction methods.The deviations vary between 1.8 cm and 15.6 cm.Overall, the lateral errors decrease with increasing water surface complexity.3. RMSE of the discrepancies between refraction corrected ALS data and terrestrial reference data in percent of the water depth.

Comparison
Overall, the effects predicted in the numerical simulation are confirmed by the experimental validation.Compared to the experimental results the simulation tends to be too optimistic.The wave induced height and planimetry coordinate errors obtained from the experimental validation are 10 times larger than the values predicted in the simulation.The main reason is the limitation of the simulations to identical forward and backward laser pulse paths here.Effects of dispersion and diffuse reflections at the water bottom are neglected.
Another reason is the difference between real and simulated water surface.The simulated wave pattern is characterized by smoother wave crests as shown in figure 4. Furthermore, the artificial character of the machine-made waves is not optimally reproduced with the oceanographic statistics based surface wave model for ocean waves, especially reflections of waves at the pool wall.
Moreover, the distribution of the water surface points differs in simulation and experiment.The simulation is based on a homogeneous distribution with 1 respectively 10 points per square meter.In contrast, the water surface representation used for the experimental validation is characterized by an inhomogeneous distribution of water surface points.Figure 8 shows an example for a typical point distribution present in the data sets.The point density varies from 0 to 8 points per square meter.The results of the experimental validation remain mean values for the whole investigation area averaging zones with high and low point density.
Assessing the results of numerical simulation and experimental validation, furthermore, the different level of detail used for the ray path modeling has to be considered.In contrast to the simulation, the refraction correction of the measurement data does not take into account the beam divergence, instead the laser pulse is treated as infinitesimal small line.

CONCLUSION
This contribution investigates the effect of wave patterns on coordinate accuracy in airborne LiDAR bathymetry.For this purpose, a numerical simulation as well as an experimental validation was carried out for a real world scenario.The comparison of the errors predicted in the simulation and the errors derived from real measurement data show a consistent tendency.However, the experimental validation reveals significantly larger errors than predicted in simulation.
Furthermore, the results indicate that the effect of wave patterns is not sufficiently considered in common refraction correction methods.Therefore, a next step could be to model local wave patterns on the basis of real LiDAR bathymetry water surface reflections.As a final result, point-wise strict coordinate correction terms can then be applied in order to increase the accuracy potential of airborne LiDAR bathymetry.

Figure 1 .
Figure 1.Refraction based on the true local water surface (blue), a assumed horizontal water surface (purple) and a assumed locally titled water surface (red) resulting in a lateral displacement dXYhz and dXYtilt and a height displacement dZhz and dZtilt.

Figure 2 .Figure 3 .
Figure 2. Empty wave pool (a) and intensity coded reference point cloud (b) acquired by terrestrial laser scanning.
The ray path modeling is realized by dividing the incident laser pulse into a large number of subbeams representing a finite footprint at the water surface.The intensity distribution within the incident laser pulse follows a Gaussian intensity profile.The refraction effects at the air/water interface are modeled by Snell's law for every individual subbeam.Our simulations are limited to identical forward and backward laser pulse paths here.Effects of diffuse reflections at the water bottom with fractions of the diffusely reflected signal accidentally being projected towards the receiver aperture are neglected.The final ground reflections are represented by the intensity-weighted centroid of all individual subbeams.In order to quantify the total effect of waves, the simulations compare laser pulse paths resulting from the refraction at the local wave-induced water surface (fig.1, blue) to paths resulting from the refraction at the horizontal (fig.1, purple) or local tilted (fig.1, red) water surface assumed in conventional correction methods.The assumption of a horizontal water surface is realized by local horizontally oriented water surface elements at different heights provided by the water surface pulse echoes.For the locally titled water surface we perform a Delaunay triangulation for all water surface points.The water surface point density is adapted to the distribution actual present in the data set.The incidence angle αtilt required by Snell's law is calculated with respect to the surface normal of the triangle intersected by the incoming laser ray.

Figure 4 .
Figure 4. Water surface profile in measurement data (a) and simulation (b).

Figure 5 .
Figure 5. Coordinate displacement error (percentage of water depth) obtained from 100 runs of the simulation tool with a flying height of 500 m, (a) horizontal water surface, (b) tilted water surface point density 3 points per square meter (p/m 2 ), (c) tilted water surface point density 10 p/m 2 .

Figure 6 .
Figure 6.TLS reference point cloud with location of the concrete base (black rectangle).

Figure 7 .
Figure 7. Deviation dZ between refraction corrected ALS point clouds (colored) and TLS point cloud (grey) at the pool bottom.

Figure 8 .
Figure 8. Water surface point density in points per square meter.

Figure 7 Figure 9 .
Figure7shows the comparison between the corrected ALS point clouds and the TLS reference point cloud at the pool bottom for one flight strip at a height of 500 m.Bottom points near the pool wall are eliminated from the evaluation to ensure that only areas with natural wave movements conforming to the oceanographic wave model are included in the analysis.Furthermore, shallow water areas are excluded, where the water depth is not sufficient for a meaningful investigation.The deviations demonstrate that the depth errors decrease with increasing complexity of the water surface representation.The main improvement results from the consideration of the local height of the water surface elements.The local surface tilt is less relevant for the depth coordinates.The differences between the ALS point cloud corrected with the simplest correction method M1 and the TLS reference point cloud clearly displays effects of the local wave pattern on the water body bottom (Fig.7 (a)).The other two correction methods leave some remaining errors as well, but less distinctive.The comparison with the density and distribution of the water surface points visualized in figure8shows that the largest deviations to the reference data occur in areas with less water surface information.The investigation of the lateral deviations between corrected ALS

Table 3
summarizes the results for all flight strips.The root mean square error (RMSE) of the lateral coordinate displacement varies in the range of 8.24 % to 11.01 %.The depth displacements are comparatively small (RMSE 1.08 % to 2.20 %), whereas the displacements decrease with increasing water surface complexity.Effects due to the different flying heights are not recognizable.