MAPPING OF WATER SURFACE LEVELS AND SLOPES WITH SINGLE PHOTON LIDAR – A CASE STUDY AT THE RIVER RHINE

Precise knowledge of water surface level heights is crucial for safe ship navigation and as basis for calibration of hydrodynamicnumerical models. While Airborne Laser Scanning (ALS) is a well established technique for topographic mapping, ALS-based water surface mapping using conventional infrared lasers suffers from the high degree of specular reflection which leads to data voids for off-nadir angles beyond 5-7 degrees. The advent of single photon sensitive ALS systems using green laser sources presents the prospect of large-area, high-resolution water surface mapping due to the high receiver sensitivity and measurement rate of such systems. Building on previous studies on subject matters, we present the results of a pilot project initiated and conducted by the German Federal Institute of Hydrology (BfG, Koblenz) at the Rhine River. Three specific test sites with varying water surface and flow velocity properties were captured on October 30 and 31, 2019 with the Leica SPL100 from flying altitudes of 3000 m, 2500 m, 1600 m, and 800 m, respectively. As anticipated, the water surface laser pulse density was high and exhibited 20-145 points/m depending on flying altitude. After quality control, strip adjustment, and point cloud analysis, three water surface classification methods were implemented based on: (i) height quantiles, (ii) point cloud segmentation, and (iii) inverse DTM filtering. All approaches featured relative and absolute water level height accuracies better than 10 cm. We conclude that Single Photon LiDAR based high resolution mapping of water surface levels and tilts is feasible when employing application specific data acquisition parameters, i.e., off-nadir angle ≤10 ◦ and flying altitude ≤3000 m.


Motivation
The availability of water surface level information is a precondition for maintenance purposes by the responsible German Waterway and Shipping Administration and for safe ship navigation within inland waterways. Furthermore, it is a valuable data source for calibrating and validating hydrodynamic numeric models, e.g., for flood simulation which is an important task for several institutions like the German Federal Institute of Hydrology (Bundesanstalt für Gewässerkunde -BfG). While water level information is sufficient for river maintenance, numerical models would greatly profit from the additional availability of precise water surface slopes.
Airborne Laser Scanning (ALS) is the state-of-the-art technique for wide-area topographic mapping which, in principle, includes acquisition of water surfaces. However, the return signal is dominated by specular reflection from the air-water-interface and high absorption in the water column. Thus, the received signals are generally weak, especially for off-nadir angles beyond 10 • . The latter especially applies to scanners employing laser wavelengths in the infrared domain of the electromagnetic spectrum (e.g., λ=1064 nm or 1550 nm), which are primarily designed for mapping topography. The penetration into the water column is minimal for theses wavelengths (Pfennigbauer, Ullrich, 2011), which is beneficial for mapping the water surface, but the dominant specular component leads to data gaps (Höfle et al., 2009), especially at very smooth surfaces like in-land lakes and reservoirs, dammed-up rivers, and river sections with laminar flow.
The Single Photon LiDAR (SPL) technology uses green laser radiation and very sensitive receivers (Degnan, 2002, Degnan, 2010, Degnan, 2016. This entails the possibility to capture weak reflections from the water surface and the topmost layer of the water column. Due to the inherent water penetration capability of green lasers radiation, water surface echoes always constitute a mixture of direct interface returns and subsurface volume backscattering (Guenther et al., 2000). While this effect is well studied for conventional Airborne Laser Bathymetry (ALB) sensors including modern high-resolution topobathymetric scanners (Mandlburger et al., 2013), only a few studies have focused on the potential of SPL-based water surface mapping and the results indicate that standard SPL flight mission parameters optimized for efficient topographic mapping (Stoker et al., 2016) lead to systematic water surface underestimation (Mandlburger, Jutzi, 2019). Based on these findings, the BfG initiated and conducted an SPL pilot project at the Rhine River specifically tailored for the purpose of large-area, high-resolution water surface mapping. In this contribution, we present the first results of this study.
The remainder of the article is structured as follows: Section 1.2 briefly summaries previous work on subject matters. Section 2 introduces the study area and data acquisition details, and describes the employed quality control procedures and initial exploratory data analysis steps. The water surface classification methods are outlined in Section 3 and the results are presented in Section 4 and critically discussed in Section 5. The main findings are summarized in Section 6.

Related Work
The general principles of single photon sensitive LiDAR are described in (Degnan, 2002, Degnan, 2010, Degnan, 2016. The technology is used for large-area topographic mapping especially in the context of national elevation programs (Sugarbaker et al., 2014, Stoker et al., 2016. The main advantage of the technology is the high areal coverage achieved by splitting a single laser beam into an array of 10x10 sub-beams (beamlets), each featuring a separate single photon sensitive receiver. The higher coverage comes at the prize of a reduced accuracy compared to conventional state-of-the-art full waveform based systems (Brown et al., 2020. While the technology is generally designed for mapping topography, the use of green laser radiation (λ=532 nm) together with conical scanning results in moderate bathymetric capabilities, which has been showcased for airborne and space-borne applications (Degnan, 2016, Parrish et al., 2019. The interaction of green laser beams with the water surface is complex and consists of mainly specular reflection from the airwater-interface itself and volume backscattering in the water column just beneath the surface (Guenther et al., 2000). Empirical studies based on topographic and topo-bathymetric laser scanners have demonstrated that systematic water surface level underestimation caused by penetration of the green laser signal into the topmost layer of the water column can reach more than 10 cm depending on the smoothness of the surface and water turbidity, which can be reduced by spatial aggregation and height quantile based statistical analysis to a few centimetres (Mandlburger et al., 2013).
In the context of SPL-based water surface mapping, previous studies have indicated that standard SPL flight configurations with a flying altitude around 4000 m and a scanner Field-of-View (FoV) of 30 • do not deliver reliable water surface estimates (Mandlburger, Jutzi, 2019). This has been demonstrated both theoretically and empirically by model calculation based on the laser-radar-equation and comparison of surveyed data, respectively. However, model calculation also revealed that more than 10 photons from water surface reflections can be expected for a flying altitude of 2000 m and a laser beam off-nadir angle of 10 • and more than 40 photons, if the flying altitude is further lowered to 1000 m.
Based on this assessment, we designed a respective data acquisition with the Leica SPL100 scanner (Leica, 2020) at the Rhine River around Koblenz, Germany, and investigated the water surface response from different flying altitudes ranging from 800 m to 3000 m employing a small, constant laser beam off-nadir angle of 10 • .

Study area and data acquisition
The BfG commissioned an SPL aerial survey, which was carried out on October 30 th and 31 th , 2019. Three specific study areas were captured from flying altitudes of 3000 m, 2500 m, 1600 m, and 800 m, respectively (cf. Figure 1). The northernmost area-of-interest (Area A) comprises the free flowing Rhine River and two impounded rivers (Mosel and Lahn), and  Table 2) features different water and land surfaces types. While the impounded river sections exhibit calm water surfaces, the surface is more riffled in free flowing part (Rhine), where ships further increase surface undulations. The dry land parts of Area A contain urban and agricultural areas, river banks, forests, and grassland. Area B exhibits two narrow river bends with expected cross sectional tilts. Area C (cf. Figure 2) is characterized by complex hydraulic conditions yielding heterogeneous river bed topography, rocks, sandbars, and narrow curves. Especially in this section, surface slopes in longitudinal and cross direction show a high variation which, in turn, requires precise modelling of the complex surface. Area C was therefore chosen as specific test site for this study.
Data acquisition was conducted with the Leica SPL100 scanner (Leica, 2020) using a 10 • scan wedge (i.e., total scanner FoV=20 • ) and a pulse repetition rate of 45-50 kHz depending on flying altitude. As each laser pulse is split into an array of 10x10 beamlets, the maximum achievable measurement rate is 5 MHz in case all channels return at least a single echo. Considering the employed constant flight velocity of 150 knots (75m/s), the pulse density (PD) solely depends on the flying  Table 1 provides an overview of the pulse density on dry land and water for the study dataset and a comparison dataset (City of Vienna) using standard SPL100 flight mission parameters described in . The table reveals that the Rhine River dataset features a higher water point density ratio of >60 % compared to 2-40 % for different water bodies of the Vienna dataset. For the Rhine dataset, the point densities generally increase with decreasing altitude. The same applies to the water-land ratio dropping from 87 % for a flying altitude of 800 m to 59 % at 3000 m. Due to the small scan angle, still a considerable amount of water points is available even from high flying altitudes. Standard SPL parameters used for the Danube River datasets, in turn, show a further density drop, which is most pronounced for the "New Danube", a relatively clear, standing water body with calm water surface. This is a first indication that the chosen application specific flight setup facilitates the task of water surface mapping.

Initial data preprocessing and data analysis
All collected data were initially processed by Leica using the HxMap Software. The trajectory was determined in postprocessing (Kalman filtering) using the GNSS corrections from the German SAPOS service (Spatial reference system: ETRS89/DREF91, realization 2016). The captured 3D point clouds were delivered in UTM projection (UTM 32N, EPSG 4647) with ellipsoidal heights.
In order the be able to directly compare the derived water surface heights with gauge levels, the ellipsoidal heights were transformed to the official German height realization DHHN2016 (EPSG:7837) using the German Combined Geoid GCG2016. At the Rhine River, an additional reference water level (GLW) is established serving as basis for nautical charts. In an additional preprocessing step, the elevations of all water points were transformed from DHHN2016 to this reference level, where the elevation values are around zero facilitating height based point data selection.

Quality control and strip adjustment
Before using the 3D SPL point clouds for water surface modelling, routine quality control procedures were applied to ensure adherence to standards w.r.t. relative strip fitting precision and absolute flight block georeferencing accuracy. Strip fitting precision check included the following steps: (i) splitting of individual strips into forward and backward looking scan lines, (i) homogenization of point density to 4 points/m 2 by selecting the the point with the highest amplitude within cells of 0.5 m edge length, (iii) calculation of Digital Elevation Models (DEM) for each (half) strip, and (iv) calculation of DEM of Difference (DoD) models for all overlapping DEMs. For the Area C block, analysis of the height deviations confirmed satisfactory results for forward and backward look of the same strip but also revealed systematic effects in the overlap area of adjacent strips. On the first glance, the height residuals indicated a small roll angle boresight misalignment.
To compensate the remaining systematic effects, rigorous strip adjustment (Glira et al., 2016) was carried out using (i) the timestamped LiDAR points in scanner coordinate system, (ii) the flight trajectory (time stamp, position and attitude), and 12 ground control patches as input. Figure 3 shows the resulting strip height differences after strip adjustment. The plot reveals that not all systematic border effects could be entirely compensated. Discussion with the scanner manufacturer confirmed that internal calibration issues (intensity dependent range correction) cause the remaining deviations. The effect is more pronounced the strip boundary due to the higher laser point density at the strip boundary along with giving preference to the high amplitude points in the above mentioned density homogenization procedure. Both the relative and absolute accuracy measures, however, are satisfactory and entirely inside the sensor's accuracy specification (relative strip fitting accuracy: 0.0±1.7 cm, absolute block fitting accuracy: 0.6±1.7 cm). The International Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences, Volume XLIII-B1-2020, 2020 XXIV ISPRS Congress (2020 edition)

METHODS
In this Section the specific methods for classifying water surface points are detailed. All methods take the potential water penetration property of green laser radiation into account (Guenther et al., 2000) and pursue different approaches to separate water surface from water column echoes. In all cases the following general processing pipeline applies: 1. Derivation of water-land-boundary (WLB) as outline of all points classified as water surface or water column (i.e., classes 9 and 64) 2. Selection of all laser points within WLB perimeter classified as water surface (class 9), water column (class 64), or unclassified (class 1) 3. Transformation of laser point elevations from ellipsoidal height system via DHHN2016 to GLW reference level as outlined in Section 2.

Statistical approach
The first method uses statistical reasoning following the general approach described in (Mandlburger et al., 2013) and consists of the following processing steps: • rasterization of laser points into 1x1 m 2 cells • calculation of height histogram per cell • classification of (potential) water surface points by selecting: the n highest points; n={5,10} the highest q % q={5,10}

Segmentation-based approach
The second approach is a refinement of the statistical method and comprises the following steps: • rasterization of laser points into cells of 5x5 m 2 • calculation height histogram for each cell • selection of 95-98 % height quantile points as water surface candidate points (i.e., seed points) • seeded region growing (search radius: 1 m, criterion: dz>-1.5 cm and dz<3.0 cm).The small negative dz limit (-1.5 cm) prohibits water column laser echoes to be classified as surface, while the larger positive limit (+3.0 cm) aims at including small wave riffle peaks in the water surface segment.
• classification of all points within segments featuring at least 75 points as water surface

Inverse DTM approach
The last approach exploits the property of decreasing point density with increasing water depth. As such, it can be seen as an inverse DTM filter problem, and water surface classification can be achieved by: • reversing the sign of the point elevations • running conventional ALS DTM filtering, e.g., based on hierarchical robust interpolation (Kraus, Pfeifer, 1998).
• re-classification of ground points as water surface

Figure 4a
shows a color coded map of the water surface model (WSM) for Area C derived from the SPL laser point cloud captured from a flying altitude of 2500 m using the statistical approach outlined in Section 3.1. The WSM generally shows spatially coherent water levels and features relative height variations w.r.t. the reference water level ranging from 0.8 m to 1.25 m. Thus, the entire height undulation amounts to 45 cm. A cross sectional water surface tilt is clearly visible in the river bend in the center of Area C. As expected, the water levels are higher at the outer bank. The tilt is pronounced for the reach around river section 551.0 km. In turn, homogeneous and horizontal water levels can be expected in the straight river reach from 549.5 km to 550.4 km.
While this is well the case for the reach between 500.2 km and 500.4 km, the section ranging from 549.9 km to 500.2 km shows a North-South-disruption, which is aligned to the borders of the corresponding flight strips. In the northern part, the same effect is visible from 551.3-551.8 km. This indicates that the employed surface estimation method is sensitive to the much higher point density at the strip boundary resulting from the circular scan pattern of the SPL100 instrument. Due to prioritizing the higher points in the surface estimation workflow, the higher point density entails higher water level estimates, which can be compensated by homogenizing the point density in an additional preprocessing step.
For the same area, Figure 4b shows the WSM derived via the inverse DTM approach outlined in Section 3.3. In general, the same patterns as in Figure 4a are discernible confirming that The International Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences, Volume XLIII-B1-2020, 2020 XXIV ISPRS Congress (2020 edition) both approaches deliver consistent results. The local water level differences between both surface estimation approaches become visible in the color coded DoD model depicted in Figure 4c and the corresponding histogram in Figure 4d. Compared to the statistical approach the inverse DTM method features a positive bias of 2.4 cm (mean/median), also indicated by the greenish color tones in Figure 4c. While the statistical approach explicitly favours higher points in a certain neighbourhood, the inverse DTM approach acts more locally and just filters water column points. For local areas with a low degree of water surface points, this leads to water surface underestimation. In addition, the strip overlap artefact clearly stands out in the height difference map, which indicates that the inverse DTM is more prone to density differences. Figure 5 depicts exemplary vertical sections of selected patches of the Mosel, Rhine, and Lahn Rivers (from top to bottom in Figure 5) captured from different flying altitudes. The red points denote laser echoes classified as water by the segmentation based approach detailed in Section 3.2. It can clearly be seen that the water point density drops with decreasing flying altitude, but nonetheless, data captured from the highest flying altitude of 3000 m still feature an aggregation of laser echoes around the water surface. The lower the flying altitude the higher the number of laser echoes from the water column. This especially applies to the 800 m dataset displayed in the lower plot of Figure 5. The gaps between the individual point columns are caused by operating the sensor outside its nominal altitude range, which leads to non-overlapping beamlet arrays. Table 2 summarizes the results of the absolute accuracy assessment for different variants of the statistical approach. For all gauging stations within the project area, the measured reference water levels in relation to GLW are compared to the SPL derived estimates and the residuals are reported for the individual water surface modelling strategies outlined in Section 3.1, namely: (i) 5% quantile, (ii) 10% quantile, (iii) 5 highest points, or (iv) 10 highest points per cell. Table 2 shows that the deviations generally exhibit a positive bias, i.e., the LiDAR derived water levels are too low. The average deviations  Table 2. Accuracy assessment of LiDAR derived water levels compared to gauge reference data; * in dammed sections, the postfixes OP/UP refer to the upstream/downstream locations of the gauge (i.e. before/after the dam)  Table 2 underline the water surface underestimation tendency ranging from 2-5 cm. The variant based on 5% of the highest points within a cell conforms best to the reference gauge water levels with median deviation of 1.4 cm and a standard deviation of ±3.4 cm calculated over all gauge stations and flying altitudes. For this variant the maximum deviation amounts to 11.9 cm at gauging station TG Lahn Lahnstein OP (flying altitude: 1600 m). While this gauge is located in a dammed section of the river, the largest deviation within the free flowing section of the Rhine River measures 5 cm at gauge station TG SanktSebastian.

DISCUSSIONS
The main purpose of the conducted pilot study was to answer the question if Single Photon LiDAR in general and Leica SPL100 point clouds in particular can be used as basis for reliable areal water surface mapping in high-resolution. Previous studies indicated that conventional bathymetric LiDAR data are better suited (Mandlburger, Jutzi, 2019) in case of using standard SPL flight mission parameters (flying altitude: 4000 m, off-nadir laser beam angle: 15 • . The current study confirmed that lowering flying altitude (800-3000 m) and scan angle (10 • ) had the anticipated positive impact for SPL based water surface detection. The captured point clouds showed a clearer answer from the water surface compared to published studies using standard acquisition parameters.
Still, as SPL uses green laser light, the inevitable penetration of the laser radiation into the topmost layers of the water column leads to slight water level underestimation of approximately 2 cm. This holds for all tested water surface classification and modelling strategies (i.e., statistical, segmentation based, and inverse DTM approach). Especially the statistical approach proved to deliver accurate results when compared to reference water levels measured at gauging stations with an overall RMSE of approximately 4 cm. Considering (i) the complexity of laser light interaction with the medium water (Tulldahl, Steinvall, 1999, Abdallah et al., 2012 and (ii) remote sensing based data acquisition from high altitude, these results are rated as satisfactory.
While all presented surface classification and modelling approaches generally performed well, the statistical approach performed best concerning reliability whereas the inverse DTM approach featured more details. The latter, however, can only be exploited if enough water surface echoes are available. One of the outstanding questions is still to reliably separate direct airwater-interface and sub-surface returns from the top most layer of the water column purely based on the captured point clouds.
Future work on subject matter will, therefore, include intensity information in the surface classification process, as higher intensities indicate returns from specular reflections from the air-water-interface. A better quantification of the water surface point classification probability could lead to spatially varying aggregation levels resulting in a water surface description with higher spatial resolution whenever enough reliable surface returns are available (wavy surface, turbid water) and to a coarser product in areas dominated by sub-surface returns.

CONCLUSIONS
In this article, we presented preliminary results of a case study aiming at areal water surface mapping in high spatial resolution based on Single Photon LiDAR. A specific flight setup was designed and carried out at the Rhine River near Koblenz, Germany, with the Leica SPL100 scanner. Building on previous work on subject matters, the main idea was to use flying altitudes of <3000 m and a small off-nadir scan angle of 10 • to increase the signal strength of weak water surface reflections and, thereby, to increase the likelihood of echo detection with single photon sensitive receivers.
The captured laser point clouds confirmed a higher level of water points compared to data acquisitions with standard SPL flight mission parameters (flying altitude: 4000 m, laser beam off-nadir angle: 15 • . Among the water points a clear density maximum around the surface could be observer. However, due to the use of green laser radiation, penetration of the signal into the water column is inevitable requiring appropriate procedures to separate water surface and water column points. We introduced three simple, yet effective water surface classification approaches: (i) a statistical, height quantile based approach using the highest 5% of the points with 1x1 m 2 raster cells, (ii) a region growing segmentation based method, and (iii) inverse DTM filtering based on hierarchical robust interpolation.
Although data processing is still ongoing, preliminary results indicate that: 1. the employed scanner settings enable SPL based water surface modelling in high resolution 2. the precision of redundantly captured water level elevations derived from overlapping flight strips or from flight strips flown in different altitudes is in the cm range 3. the inevitable water surface underestimation due to partial penetration of green laser radiation into the topmost layer of the water column is less for flight strips flown in lower altitudes (e.g., 800 m) 4. longitudinal water surface slopes can be captured reliably and cross sectional tilts in narrow river bends are detectable in the point clouds 5. the mean absolute deviations compared to independent reference gauge levels are <5 cm Based on this preliminary assessment, the Leica SPL100 shows high potential for enabling high-resolution water surface modelling when application specific flight mission parameters are applied. Precise water surface data may be used as basic geodata for nautical charts, reference data hydrodynamic-numerical models, and the like. Data processing of the Rhine River pilot project is still work-in-progress. Future investigations will focus on a better understanding of the optimum flight mission parameters balancing water surface precision and flight mission efficiency.