COMPARISON OF SUBAQUATIC DIGITAL ELEVATION MODELS FROM AIRBORNE LASER SCANNING AND IMAGERY

The paper describes and compares the workflows and results of generating digital elevation models (DEMs) of underwater areas from airborne laser scanning and aerial stereo images. Based on a combined laser scanning/image data set of an artificial lake, both methods are described and pros/cons are highlighted. The authors focus on the final results, especially on accuracy, completeness and spatial resolution of the underwater DEM’s. Further, practical aspects of processing and complexity of both methods are highlighted too.


MOTIVATION
Lasers scanning and photogrammetry are two well established methods for Digital Elevation Model (DEM) acquisition in all scales. Before the advent of laser scanning systems, stereo photogrammetry was the method of choice for this task. Over the past decades, laser scanning was established as an effective tool for DEM acquisition. Especially the multi-target technology as well as full-wave-form analysis  tipped the balance in favour of laser scanning. Nevertheless, the introduction of highly automated techniques for image orientation (e.g. Structure from Motion, SfM) and dense pixelwise matching led to a renaissance of photogrammetry in the last decade. As shown in several projects (Mandlburger et al., 2017), both methods can be applied complementary in order to improve the final results.
Airborne laser bathymetry established its stand as a standard method for capturing submerged topography of shallow water Figure 1. Study area Autobahnsee-a manmade freshwater lake near Augsburg, Germany. * Corresponding author bodies, like riverbanks or shore lines (Wozencraft, Lillycrop, 2003, Kinzel et al., 2013, Song et al., 2015. For the same task, photogrammetry can be used in the classical way. Both methods have to take the refraction into account while measuring through refracting surfaces -in this case water. Crucial for modelling the refraction is precise knowledge of the water level height and orientation as well as the refractive index. The refractive index for a specific water body can be estimated from salinity and temperature (Quan, Fry, 1995). More demanding is the determination of the shape of the water surface, especially in wavy conditions. Further, effects of water flow, causing local sinks and stagnations, have to be taken into account. Nevertheless, in the case of calm water, the water surface can be approximated well as a horizontal plane. Therefore, only the water level height has to be determined. This can be done via local gauge measurement or simultaneously within data processing.

DATA SET
Laser and image data were acquired simultaneously on 9 April 2018 during a flight over a lake in the flood plain of the river Lech near Augsburg, Bavaria, Germany (see Figure 1). The main substrate on the lake ground is coarse grained gravel. The airplane was equipped with a topo-bathymetrical laser scanner (Riegl VQ-880-G) together with two 100 Mpix cameras (IGI DigiCAM 100). One was a standard RGB-camera, while the second one was operated in a monochromatic mode together with a Coastal-Blue filter (λ=400-460 nm) . Altogether four stripes were acquired at flight speed of 100 knots (50 m/s).

Image Data
The image block consists of 125 images captured from two different flying altitudes; 64 images from 610m resulting in a ground sampe distance (GSD) of 5.6 cm and 61 images from 450 m with a GSD of 4.2 cm. For each height, two strips were acquired with an overlap along track of 90% and cross track of 60% (see Figure 2). Therefore, the redundancy is quite high, e.g. the centre of the area was captured in 40 images. As mentioned before, the used camera was a IGI DigiCAM 100 medium format camera with a resolution of 10608x8708 Pixel (100 Mpx) together with a 50 mm lens. The camera itself was developed from a PhaseOne iXU-RS 1000 . For the tests, the RGB images were processed exclusively in favour of the Coastal-Blue data, because of the more suitable characteristics for automatic image analysis . For geo-referencing GNSS and INS data were acquired during the flight. Additionally, 10 ground control point targets were installed and measured via RTK-GNSS (see Figure 2).

Image Orientation
In a first step, the whole block was processed in Pix4D Mapper (www.pix4d.com) with geo-referencing based on control points. In the same software, a point cloud together with an orthophoto mosaic was generated. While the DEMs can be assumed correct for land areas, the submerged surface points were falsified by refraction. Generally, water depths in these areas are underestimated when refraction is not taken into account (Maas, 2015). So, for a correct calculation of point heights a compensation of refraction is essential. Thanks to previous works in this field, a dedicated software was available. This multi-media bundle adjustment software was developed at the Institute of Photogrammetry and Remote Sensing, Technische Universität Dresden (Mulsow, 2010).
Processing of the whole block in MatchAT (Trimble/ Inpho) provided initial values for interior and exterior orientation as well as image point coordinates of control points (manual) and tie points (automated), respectively. For classification of submerged and land points, the heights of the 3D tie point coordinates were usedwith the water level height as discriminator. From manual image measurements of the shore line, a water level of 509.5 m could be determined. All points below this level could be classified as submerged. Altogether, ~17.000 image measurements of ~1000 tie points were provided for the multimedia bundle adjustment, of which ~1900 image coordinates related to ~120 submerged tie points.
Experiences on similar projects had shown that a reliable determination of the water level inside the bundle process is not possible from block of nadir images (Mulsow, 2018). Therefore, the lake level was defined as a horizontal plane with a fixed height of 509.5 m. A reliable determination of refractive surfaces inside a bundle adjustment requires convergent images as well as relatively low impact angles of image rays on to the refractive surface. As shown in  for close range image blocks with proper block geometry, the surface parameters as well as refractive indices can be determined with high accuracy.
Bundle block adjustment was done in two parameter configurations with different treatment of underwater points: I. Simultaneous determination of all unknowns (interior and exterior orientations), underwater points treated as single-media points (without modelling of refraction)

II.
Simultaneous determination of all unknowns (interior and exterior orientations), underwater points treated as multi-media points (considering refraction) Both adjustments were processed without gross error detection, which was done in advance. As mentioned before, both runs were done with the same input datajust for the underwater points the refraction was either neglected (configuration Iconventional bundle adjustment) or modelled (configuration IImultimedia bundle adjustment).
From Table 1 it's obvious that the modelling of refraction has virtually no impact on the accuracy values (back projection error s0, RMS of residuals of image measurements, RMS of the standard deviations of the estimated object point coordinates). It is noted that the RMS values reported in Table 1 are calculated from the (a posteriori) residuals rather than from the covariance matrix of the bundle block adjustment. The positions of tie points are close together for both adjustments (RMS 0.1 m in x/y, note noted in Table 1). As expected, the heights differ significantly for submerged areas (see Figure 3).  Anyway, the exterior orientation parameters as well as camera parameters are similar for both runs. The differences are smaller than their accuracy, which means that a conventional bundle adjustment is suitable for block orientation in this example. It can be assumed, that vertical oriented image blocks showing shallow water bodies can be processed as standard blocks.

Parameter
The reason lies in the homogeneity of refraction effects: Therefore, image rays of submerged points are refracted in a similar way. Further, the water depth to flight height ratio is quite small, which means that only a small part of the image ray is refracted. In the end, both effects result in a relatively small distance of (refracted) image rays from their adjusted (submerged) intersection point. According to this, standard software can be used for orientation of such image blocks Anyway, no definitive recommendation can be given without further investigations. For DEM extraction a dedicated intersection routine with refraction compensation should be used (Mandlburger, 2018). As shown in (Mulsow, 2010) such a routine can be easily implemented, in contrast to a resection-or bundle-adjustment software.
Table 1 (Setup II) shows a drop of accuracy of just factor 2 (RMS from adjustment) for land and underwater points. This can be seen as far too optimistic, because for calculation of inner accuracy values the refraction index as well as the water surface parameters were treated as fixed. Nevertheless, small local variations of these parameters due to wind and temperature can be assumed. A comparison with check points would give a more realistic estimation of height accuracy. Unfortunately, no such reference data were taken in the field. From similar projects (Mulsow, 2018), a drop of accuracy of underwater points of factor 4 (against land points) can be assumed. Therefore, a height accuracy of 10 cm instead of 5 cm can be estimated.

DEM from image data
The DEM extraction follows the common procedure: -image matching in stereo pairs -computing of object coordinates via forward intersection. First, image pairs were defined automatically from the low altitude image subset. In order to achieve large intersection angles, pairs with an overlap of ~60% were chosen. Due to high overlap along track, a number of neighbouring images could be neglected. The average distance of stereo partners was about 5 images. An overlap of neighbouring image pairs of 60% was defined too. In the end, only half of the low altitude images were used.
Stereo image pairs were rectified with normal geometry in order to minimize y-parallaxes. In contrast to single-medium case, the epipolar lines are not straight for underwater points due to refraction. However, thanks to the vertical imaging direction as well as the relatively shallow water (compared to the flying altitude), this effect is negligible for image point matching.
The matching was implemented as a pyramid search approach, starting from a reduced resolution of factor 5 up to the full resolution. For each step, the image was divided into 75x50 raster cells. For each cell, the best feature point (Harris-operator) was found (Harris, Stephens, 1988). The Harris points were localised in the stereo partner image and the subpixel position was measured via Least Squares Matching (Gruen, 1985) with a patch size of 21x21 pixel, x/y shift and scale parameter. From parallax differences, a disparity map was created for the actual resolution step. The disparity map was resampled for the next resolution step and gave initial value for the following matching process. After reaching full resolution, a final matching was done with full parameter set in order to achieve best results and to overcome refraction effects. Based on matching results, all object points coordinates were computed via conventional forward intersection.  In order to identify underwater points, the points were analysed by their heightspoints below the pre-defined water level were considered as submerged points and vice versa. Finally, the underwater point coordinates were computed via multimediaforward intersection (Mulsow, 2018).

Analysis of results
As mentioned before, the focus of this study laid on the DEM quality of submerged areas. As shown in Figure 4, underwater feature points were found mainly on the shore and on border zones of vegetation. This illustrates the main drawback of image based DEM generationthe dependency on (image) texture. Therefore, homogeneous areas show either no or false measurements. Nevertheless, for this project points up to 4 m below water surface (max. depth 4.5 m) could be determined with good reliability (see also Figure 5).
The precision of point coordinates is not homogenous for the whole area. As expected, the imaging quality drops with increasing water depth. Therefore, the precision drops with increasing depth too. Without reference data this effect can only be estimated.
The high overlap inside the image block (along-and cross-strip) results in a high overlap of DEM's from different image pairs. Therefore, the inner precision can be estimated by comparing the height data from overlapping DEMs. For this purpose, two DEM's from different strips were transferred to cubic-splinesurfaces and rasterized to 0.5 m grids. The comparison of node heights shows a good fit of both models (RMS 4.4 cm, average distance in mm-level, no offset). Nevertheless, areas with only a few feature points (low image contrast) show significant differences (max. 22 cm). However, the precision of underwater points given in Table 1 could be confirmed in general.
In theory, the geometric conditions for forward intersection is more stable for land-than for underwater points due to refraction. The image ray intersection angle in the denser medium (water) becomes smaller (Maas, 2015), thus degrading the accuracy.
For a visual test, the surface model was intersected with the water surface and the resulting line was projected was plotted in the orthophoto. As seen in Figure 9, the calculated shore line follows the real one quite closely. This especially holds true for open areas as depicted in the detail in Figure 9, whereas larger deviations can be observed in vegetated areas.

Orientation of laser scanning data
Together with the image block, an additional data set was acquired simultaneously by a topo-bathymetric LiDAR system (Light Detection And Ranging). The whole dataset included the trajectories of 4 strips (GNSS/INS), sensor data (range, angle, calibrated amplitude, reflectivity per echo from online waveform analysis (Pfennigbauer et. al., 2014), and several other per-point attributes (echo number, scan angle, etc.). For registration, a rigorous strip adjustment including time dependent trajectory correction Glira et al., 2016) was performed. Quality control was based on the residues between the strip heights in smooth areas. The robustly estimated standard deviation (σMAD) of the remaining height differences between overlapping flight strips amounted to 1.0 cm and the residuals showed no bias after adjustment.
Absolute orientation of the laser flight block cloud was established via rigid body transformation based on the photogrammetrically derived point cloud by means of the Iterative Closest Point (ICP) approach (Besl andMcKay, 1992, Glira et al., 2015). Apart from negligible rotation components, the main part of the laser-to-image transformation was covered by a 3d-offset (dx=-0.087 m, dy=-0.028, dz=+0.514). The standard deviation of the unbiased residual height deviations (mean dz=0.001 m) between the image block and the transformed laser and block measured 0.036 m.

Determination of water surface and refraction correction
The laser ray is refracted while travelling through the water surface. Further, the propagation speed is reduced in water. Both effects are formulated in the Snell's law and depend on the refractive indices of air (nA~1.0) and water (nW~1.33). Therefore, the uncorrected laser points appear too deep for water bodies.
For correction of time of flight and beam refraction, a digital model of the water surface is essential. The main advantage of laser bathymetry is the ability to receive reflected signals from the water surface as well as from the ground below. Because of the high amount of specular reflection at the water surface and the general water penetration capability of green laser radiation, the backscattered laser signal is a mix of surface reflection and volume scattering in the first centimetres of water column (Guenther et al., 2000).
For determination of the water surface, the approach by Mandlburger et al. (2013) was chosen. Based on a manually defined initial height, the laser echoes were statistically analysed in 10 m by 10 m cells. For each cell, the representative water level was calculated as the 99% quantile of all height values within each cell. On the one side this robustly eliminates potential above surface outlier points while at the same time minimizing the water level underestimation due to the volume backscattering. The water surface model showed maximum undulation of 15 cm with systematically higher elevations around the island in the southern part of the lake. The inclusion of occasional land points in the transition zone between water and land led to a slight water level overestimation in the shore areas. Therefore, the maximum water level height was limited to h=509.10 m, while the mean height was 509.00 m. The mean laser derived water level height is 5 cm below the value from image data.
With the water surface DEM on hand, each echo was corrected by intersecting the laser ray with the water surface and by calculating and applying the refraction and run-time correction for the part of the laser beam in water. Therefore, each raw underwater point coordinate saw corrections in position and height. The whole correction process was carried out in the software system OPALS (Pfeifer et al., 2014).

Filtering and DEM
After pre-processing, the point cloud was filtered for terrain points (on land and submerged) no-terrain (low, mid, high vegetation, buildings etc.) via hierarchical robust filtering . An additional classification was carried out to separate lake floor, water surface and points in the water column based on the water surface model. All terrain points (on land and submerged) were used as an input for a regular DEM with 50 cm spacing. While the laser point density of 25 points/m 2 would allow a lower spacing, the effective spatial resolution is limited by the size of the laser footprint of 50-60 cm Beside the quality of the reconstructed flight trajectory (GNSS/INS) and the accuracy of the scanning device, the uncertainty of water surface orientation and height has to be taken into account for estimation of accuracy potential of the derived point cloud. From similar projects, a height accuracy of at least 10 cm can be assumed. Due to redundant coverage with data from 4 overlapping flight strips, the relative DEM precision can expected to better by a factor of 2 (5 cm).

COMPARISON OF RESULTS
The comparison of DEM's from both sources showed a high degree of similarity for vegetation free shore areas as well as for zones with sufficient texture (see Figure 8). Therefore, discrepancies larger than 20 cm were excluded in the analysis. The overall coverage of laser DEM is significantly higher than image based DEM. From this, the advantage of active method (laser) against passive image based when capturing low-texture areas becomes obvious.
Due to the lack of independent reference data, the absolute accuracy of point coordinates below water level could not be assessed However, comparison of the shore lines independently derived from both DEMs and water surface models allow a visual estimation of accuracy. The shore line derived from the laser DEM follows the real one (see Figure 9) in the range of a few decimetres. It has to be kept in mind, that the gradient of the shore zone is quite small and therefore the same applies for the intersection angle with the water surface. The same can be observed for the image DEM for vegetation free shore zones. By contrast, shore line sections hidden under vegetation area shifted in direction of the lake. This highlights the main advantage of laser based DEM extractionthe multi-target technology in general and as full waveform analyses in particular allow a look under vegetation. Figure 9. Shore lines calculated from intersection of water surface with image-DEM (blue) and laser-DEM (red).

CONCLUSION AND OUTLOOK
Both methodslaser bathymetry and two-media photogrammetryare competitive tools for acquisition of underwater topography. Nevertheless, laser scanning as well as imagery have their pros and cons regarding complexity, accuracy and completeness of the final result. As for all optical measurement techniques, the optical characteristics of passed media, in this case air and, above all, water, are crucial for measurement quality. Both methods are strongly influenced by dispersion and turbidity. When it comes to accuracy, a clear overall recommendation can't be given, because of the similar inner accuracy (5 cm vs. 4.4 cm, comparison of overlapping strips and models). As mentioned before, laser scanning is independent from texture like photogrammetry. Further, the full waveform analysis allows a look under vegetation to a certain degree. On the other hand, the technical effort as well as complexity is far lower when using cameras for data capture.
During image data processing, some possible improvements could be identified. When planning such measurement campaigns, the requirements of two-media photogrammetry should be taken into account. In order to achieve larger intersection angles, optics with large opening angles should be used. Oblique imagery could further improve the intersection geometry. However, the camera axis should not deviate from the direction vertical too far, otherwise total reflection could occur and wave effects could degrade the imaging quality. For reliable geo-referencing and thorough verification of final data, control points should be marked under water. This recommendation is valid for both laser scanning and multimedia photogrammetry. Further research should be invested for automatic extraction of shore line from images. From differences in brightness and colours (see Figure 9) an automatic identification of the borderline between dry and wet areas are feasible (Kroehnert et al., 2017 and. In case of calm water bodies, the shore line should be horizontal. This fact could be used as a constraint in the bundle block adjustment. Highest degrees of automation are to be expected using multispectral data including infrared channels due to the high absorption rate of water. This successfully demonstrated both images (Sivagami and Jayanthi, 2019) and laser scans (Morsy et al., 2018).