# Full-waveform LIDAR fast analysis of a moderately turbid bay in Western France

Patrick Launeau, Manuel Giraud, Marc Robin, Agnès Baltzer

► To cite this version:

Patrick Launeau, Manuel Giraud, Marc Robin, Agnès Baltzer. Full-waveform LIDAR fast analysis of a moderately turbid bay in Western France. Remote Sensing, 2019, 11 (2), 10.3390/rs11020117. halshs-02319241

HAL Id: halshs-02319241

https://shs.hal.science/halshs-02319241

Submitted on 10 Nov 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

Article

# Full-Waveform LiDAR Fast Analysis of a Moderately Turbid Bay in Western France

Patrick Launeau<sup>1,*</sup> (ID), Manuel Giraud<sup>1</sup>, Marc Robin<sup>2</sup> and Agnès Baltzer<sup>2 ID</sup>

Received: 25 November 2018; Accepted: 6 January 2019; Published: 10 January 2019

<figure><img src="image_3.png" /></figure>

check for
updates

**Abstract:** In shoreline monitoring, only topo-bathymetric light detection and ranging (LiDAR) can map large corridors from aerial dunes to sandbanks in shallow water. Increasing turbidity masking the formation of 532 nm laser beam echoes on the sea bed makes this challenging. Full-waveform recording all the laser beam damping functions, a turbid water column can be seen as an accumulation of layers forming a single continuum and a distinction can be made between signals ending at the bottom down to a depth of 10 m. In practice full-waveforms are converted by laser beam tracing an image cube with a grid of 1-m-wide pixels and a 0.15 m range resolution storing the mean intensities returned along incident angle. The first derivative of a wide Gaussian filter serves to delineate the full-waveform range limits. Because of turbidity current heterogeneity and complexity of multiple layers radiative transfer model, a drastic simplification is applied by normalizing the cumulative full-waveform to 1, transforming each pixel of the water column into spectrum of intensity ranging from 0 to 1 from bottom to top. A transposition between range and bottom-top information facilitates the water index correction below the sea water level provided by 1064-nm discrete echoes. All echoes remain accessible with maximums of cumulative full-waveform third derivative.

**Keywords:** full-waveform; airborne LiDAR; bathymetry; shoreline monitoring; turbidity

## 1. Introduction

Coastal areas are undergoing many geomorphological changes under the drivers of sea-level rise, negative sediment budget, incident wave energy and a variety of human interventions along the shoreline. There is a need to document the dimensions and rates of change to inform coastal management and to help define options for responding to the evolving environment.

Numerous studies have assessed shoreline change [1–5]. To do this, many shoreline proxies have been used [6]. For example, proxies such as dune vegetation limit or dune toe limit for sandy coasts, cliff base or seaward vegetation limit of the cliff top for rocky coasts, and the base of the coastal defenses protection structure can be used to monitor shoreline dynamics [7,8].

One can make a precise assessment of coastline mobility using such proxies. Unfortunately, this is not always enough to fully understand the observed dynamics. In fact, the simultaneous consideration of the mobility of foreshore and nearshore sediments in shallow coastal water, which would allow scrutiny of the entire sediment cell, is missing.

The observation of the nearshore and foreshore bottom is costly and time consuming. A recent alternative is the use of airborne LiDAR offering high performance when compared to multibeam echosounders. This resulted in the beginning of a national survey project on mapping the French coastal areas (sea and land) [9].

<footnote n="1">¹ Laboratoire de Planétologie et Géodynamique, UMR CNRS 6112, Université de Nantes, 2 rue de la Houssinière, 44322 Nantes, France; manuel.giraud@univ-nantes.fr</footnote>

<footnote n="2">² Littoral Environnement Télédétection Géomatique, UMR CNRS 6554, Université de Nantes, Campus Tertre, 44312 Nantes, France; marc.robin@univ-nantes.fr (M.R.); agnes.baltzer@univ-nantes.fr (A.B.)</footnote>

<footnote>* Correspondence: patrick.launeau@univ-nantes.fr; Tel.: +33-2-51125266</footnote>
In order to conduct cost effective, time-efficient monitoring of the bottom of shallow and turbid
coastal waters, there is a need for adapted methodologies to improve remote detectability through the
water column.

The development of airborne LiDAR bathymetry is summarized by Guenther et al. [10] in a presentation of the Scanning Hydrographic Operational Airborne Lidar Survey (SHOALS) system [11,12]. Almost 20 years later we are still using the Titan system from the same company, Teledyne Optech, capable of recording the full-waveform (FWF) of the laser beam return of a water surface with an infrared (1064 nm) laser beam and from the bottom of shallow water with a green (532 nm) beam resulting from a wavelength splitting in two of the same laser source. Most studies like [13] characterize the amplitude of the 532-nm full-waveform of shallow water by a first return at the air-water interface and a final lower, but rarely strong, return at the sea bottom. These are proportional to rock, sediment and vegetation cover reflection with an intermediate water volume backscattering proportional to various conditions of water color studies [14]. These conditions are (i) the water column depth and its concentrations in colored dissolved organic matter (CDOM); (ii) nonalgal particles and phytoplankton modifying both absorption and diffusion coefficients; and (iii) the interface geometry like sea waves at the upper interface and ripple marks at the bottom of the sea. In this radiative transfer model approach of the bathymetry the bottom is individualized as a specific echo whose calibration of the amplitude eventually allows the detection of the bottom type for [13]. In such cases bottom reflectance corrections are required which involve bidirectional reflection distribution function (BRDF) [15] considerations. These not only include bottom scattering but also water surface specular hot spot features unrelated to water composition and quite frequently due to the transmission and detection optics located coaxially as stated by [13]. As observed in many LiDAR studies [13,16,17], a steep slope can produce pulse stretching proportional to surface orientation depending on the footprint as studied by [18] with FWF Gaussian-like decomposition [19,20]. Although it has a minor effect [10] wave can also matter and recently Westfeld et al. [21] proposed a correction of coordinate errors induced by ocean wave pattern blurring the bottom images.

In shorelines with spatially heterogeneous daily sediment remobilization during tide flows, the increasing turbidity tends to merge intermediate and bottom echoes in a single continuous FWF function. At the water surface, wind corridors due to local coastal relief can also generate a strong local hot spot giving the greatest intensities on water interfaces. The calibration of water and bottom sediment optical properties being particularly difficult in such heterogeneous environments, we looked for an alternative approach to bathymetry without intensity calibration. The proposed method considers only the length of light interaction with a medium, whatever its composition. In this model, we only measure distances between the beginning and the end range of a signal which can be represented by a cumulative curve normalized between 0 and 1. By so doing, dark (absorbent) and bright (reflective) medium would give the same function with a distance of between 0 and 1 varying only with the range of diffusion properties. This was tested by [22] to map the low vegetation range where top and bottom did not allow the individualization of separate echoes. Because of signal stretching on the gradient this method works at the smallest possible pixel size without neighborhood contribution. In such a model any grass tuft is a tiny reflector and the resulting LiDAR FWF is a collection of an infinite number of tiny returns stacked on top of each other forming a single echo with upwards or downwards damping function which should not be confused with the laser pulse emission asymmetric waveform. Similarly, heterogeneous layered turbidity currents can produce a multitude of intermediate water interfaces. Thus, this work explored the variations of diffusion and backscatter of a laser beam along the water column, in the direction of the angle of incidence, from the surface downwards to the end of the signal on or near the bottom. The water surface can therefore be dark or bright and the bottom may correspond to a continuous progressive attenuation of FWF.
## 2. Materials and Methods

### 2.1. Study Area and Data Acquisition

The bay of La Baule is located on the Atlantic coast of France near Saint-Nazaire, northwest of the Loire estuary (Figure 1a). Its bottom gently slopes (less than 0.2%) seaward to a maximum depth of –6.8 m CD (Chart Datum see Table 1) in the mid-southern section. It consists mainly of bottom fine sediment with a few rocky flats to the west. The sediment mobility is quite active with sandy accretion in a westerly direction. There is a large nearshore with a sand bank in the central part of the bay which is an east–west oriented longshore swash bar. There are also a few nearshore narrow spit-like sand banks oriented to the southeast in the eastern part of the bay. This bay is located just north of the mouth of the Loire estuary, which partly explains the medium (30 mg L⁻¹) to high (maximum of 70 mg L⁻¹) turbidity level of the bay water, depending on tidal currents, flow of the Loire river and prevailing winds.

<figure><img src="image_4.png" /><figcaption>**Figure 1.** Location of the study areas: (a) discrete 1064 nm first echoes DSM, (b) synchronous orthophoto; (c) minimum, maximum and mean intensities of 1064 nm echoes; (d) minimum, maximum and mean intensities of 532 nm echoes; (e) discrete 532 nm first echoes DSM and (f) discrete 532 nm last echoes DSM (with water index correction provided by Lidar Mapping Suite (LMS) [23]); red +1 and +6 are sample locations of data used for presentation of the method.</figcaption></figure>

The Teledyne Optech Titan dual wave LiDAR managed by Nantes Rennes university platform can record a FWF on one selected channel at a time. It was installed to record the 532-nm signal
for bathymetric purposes. The airborne campaign was conducted at an altitude of 400 m with the cooperation of GEOFIT and PIXAIR companies over two days, 6–7 October 2017. The scan angle of the rotating mirror was ±14° giving a field of view (FOV) of 28° and other details of data acquisition are given in Appendix A.

A 1 m × 1 m per pixel Digital Surface Model (DSM) of the first 1064 nm returns envelope was generated for a quick overview showing tide variation between swath and days (Figure 1a, see also high resolution figures in Supplementary Materials). Along with the Titan system and a JAI Integrated 5MP camera, GEOFIT Expert Company, 1 route de Gachet - CS 90711 – 44307 Nantes France, also provided a 1 m × 1 m pixel orthophoto of the area showing a characteristic brown colored water with turbidity variations. The 532-nm beam penetrating the water light path and speed are modified at the air–water interface, which is corrected by LMS [23] in final 532-nm echo files, not in FWF. Figure 1c,d present the minimum, maximum and mean of pixel size discrete intensities of 1064 nm and 532 nm returns respectively, water index correction is included for 532 nm data. All diffusing materials, vegetation thickness of 1064 nm data in Figure 1c and water column of 532 nm data in Figure 1d, appear with a false green color because of a larger scattering of their returns. By contrast, compact materials appear magenta. An incorrect setting explains the particularly strong hot spot effect along the nadir of the 1064 nm strip lines, Figure 1c, on all wet areas during the first airborne day. Wet sand with a tiny water film on the top produced the strongest returns saturating the LiDAR recording. Most bathymetric LiDAR on the market use a tilt of 6° between the 532 nm and 1064 nm beam to facilitate the penetration of the 532-nm laser beam in the water and its reflection on the top of the water with 1064 nm, but the absolute orientation of both laser beams with the ground, driven by the aircraft attitude, is also important, as experienced with the excessive 1064-nm hotspot of Figure 1c. Fortunately, this defect had no impact on the 532-nm data, Figure 1d. The DSM of the first 532 nm discrete echoes does not always catch the water surface in Figure 1e. This highlights the importance of water level detection with discrete 1064 nm echoes (Figure 1a). The DSM of the last 532 nm echoes used to give the bathymetric data, including water correction provided by LMS in Figure 1f, only accessed the most translucent part of the bay. This is mainly on its northern part, sheltered from the prevailing west to east wind parallel to the limit of bottom detection (Figure 1f). This is also correlated with a more visible perpendicular wave pattern on the second day of data acquisition, in the southern area where bottom detection is more difficult.

## 2.2. Methodology

### 2.2.1. Full-Waveform Data Registration and Ray Tracing

Like many authors who have studied the diffusion of light through the thickness of vegetation [24–29], a 3D (x, y, z) voxel grid is used to resample all FWF in 401 channels with a vertical range resolution of 0.15 m acquired every nanosecond on a 2D pixel resolution of 1 m × 1 m. Nevertheless, in order to preserve the same light path in each pixel as discussed by [22], a ray tracing procedure is used to build a pseudo hyperspectral 3D cube image on an image plane at elevation 0 m with the same x and y pixel dimensions using a R LiDAR range dimension, instead of hyperspectral wavelength, corresponding to the laser beam light path until it hits the ground or sea bottom with a tight incident angle. Following [30] in a full-waveform LAS file, each discrete point is associated with a three-dimensional vector and a table of full-wave points. The vector gives the direction of the reflected laser pulse. A fixed range grid to find the minimum and maximum ranges was used. The resolution is given by the temporal spacing of waveform points retrieved from the LAS file. For each discrete point, we used the associated three-dimensional vector to compute a simple line–plane intersection with a plane at a zero-meter elevation. At the projected point, on the (x, y) grid of the image plane, we first accumulated the contribution of each waveform incident angle in the image of incident angle. Then, on the same (x, y) grid, we accumulated along the waveform vector subdivided into 401 channels of 0.15 m range, the contribution of the associated waveform at each range into the selected range channel.
Finally, we averaged channel by channel all contributions in each pixel. From a flight line altitude of 400 m, a laser beam hit 1-m-wide pixel at 0 m elevation with a local FOV of 0.3°. Consequently, all FWF averaged per 1-m-wide pixel correspond to an angular cone of 0.3°.

In mosaicking, the images are tiled on top of each other, without averaging overlapping images, for the preservation of the tight scan incident angle attached to each pixel. Alternatively, it is also possible to favor the pixels having the smallest incident angle in mosaic to minimize image border effects. In such cube images, the R dimension stores the variation of intensities along the incident angle of a vector length R and a conversion from (x, y, R) image cube to a 3D (x, y, z) voxel grid requires a resampling, redistributing the intensities between pixels and channels with respect to their incident angles. This geometry preserves the light path for a perfect matching with hyperspectral data as discussed in [22].

### 2.2.2. Full Waveform Data Normalization

This procedure was defined by [22]. Instead of calibrating the waveform amplitude, as discussed in the introduction, we defined a relative intensity independent of the amplitude to measure only the scattering of the FWF without any a priori knowledge of its shape or composition of the interface at the origin of a return. Then, a FWF can be made of an infinite number of functions distributing relative intensities anywhere in the range of the FWF. Instead of focusing the analysis on upward diffusion like [22] we focused on downward diffusion for bathymetric purposes thus maintaining the information of absolute altitude. The key issue becomes the detection of the beginning of the meaningful FWF signal which is a limit between the continuous FWF downward damping function and the blank noise of a baseline. Six sites located in Figure 1f are used to illustrate the data processing step-by-step. As shown in Figure 2a, with the example of site #1 (Figure 1f), the FWF are recorded with an offset of relative intensity 200 DN (with DN for digital number unit). But in practice this baseline can rise above 200 DN as in the example in Figure 2a or drop down to 180 DN in noisy spots. It is therefore necessary to determine the effective baseline level for each pixel with a longer record of signal down to -20 m (in French geodetic system NGF for Nivellement Général de la France) out of the bathymetric range. An example of asymmetric laser pulse emission, not to scale, is also presented in Figure 2b for comparison. The following processes use successive combinations of long convolution low-passes and short gradient high-passes as designed by [31]. It uses old mathematical morphology principles [32] of the intercept method adapted to edge detection of minerals in rocks in grey levels [33] and transposed in this work on FWF dimensions looking for turbid water layer imbrications instead of sandstone bedding.

First, FWF beginnings (at range $R_b$) and endings (at range $R_e$) are defined, Figure 2a, by the identification of the first transition from 0 to a positive value and the last transition from a positive value to 0 to prevent false detections and artifacts at both ends of the FWF. For instance, the low-pass {0.11, 0.22, 0.34, 0.22, 0.11} applied to the FWF still displays sharp bounds (1) in Figure 2a. Next, a first gradient {-1, 0, +1} followed by the double application of the {0.11, 0.22, 0.34, 0.22, 0.11} low-pass (equivalent to {0.0121, 0.0484, 0.1232, 0.198, 0.2366, 0.198, 0.1232, 0.0484, 0.0121}) provides the first smoothed derivative of the FWF (dFWF) and the second one {+1, 0, -1} also completed by the same two {0.11, 0.22, 0.34, 0.22, 0.11} low-passes applied to dFWF, provides the second smoothed derivative ddFWF, whose local maximum is greater than the threshold chosen by the operator, provides a selection of echoes. A 3-m full-width at half-maximum (FWHM) Gaussian low-pass filter modified from [34] using j continuous weighting coefficients (up to 60 in this work) sliding on the center of any FWF range channel $R_i$ of Equation (1) is applied in the FWF interval $R_b$ to $R_e$ (Figure 2a):

$$ \text{FWF}(R_i) = 0.5 \cdot \sum_{R_j=R_b}^{R_e} \text{FWF}(R_j) \cdot \exp\left(\left(\frac{R_j - R_i}{\text{FWHM}}\right)^2\right) \quad (1) $$
Its first derivative {−1, 0, +1} is then used to determine a shorter segment of meaningful signal right above the baseline estimation. In practice a forward scan of the first derivative of the 3-m FWHM Gaussian filtered FWF, to the first value above the $DN_b$ threshold chosen by the operator, gives the beginning $R'_b$ and a backward scan up to the first value below the same but negative threshold $DN_e$ gives the end $R'_e$ of the shorter meaningful FWF (Figure 2c).

These new bounds ($R'_b, R'_e$) are used to build the cumulative FWF (CFWF) focusing on the useful signal. Its division by its last value gives the normalized CFWF (NCFWF) smoothed with one {0.11, 0.22, 0.34, 0.22, 0.11} low-pass. As shown by [22] it is also possible to calculate the normalized centered CFWF (NCCFWF) by using the height of the maximum values of the FWF as a new common reference, aligning all main echoes on the same level zero, facilitating the distinction between positive upwards and negative downwards damping function of the FWF by signal diffusion or backscattering, whatever their altitude or bathymetry. The smoothed NCFWF varies from 0 at the bottom of the water column to 1 at the top, with height varying on the abscissa which is not convenient for bathymetric correction of a FWF image cube. Therefore, we define a transposition of NCFWF in range versus diffusion proportion from 0 to 1 forming a new function $NCFWF^T$ made of 200 channels ranging from 0.05 to 1 in Figure 2d. The empty gap between transposed values of $NCFWF^T$ is simply filled by the closest positive values weighted by their distance to each empty channel. A low-pass {0.11, 0.22, 0.34, 0.22, 0.11} smooths the final $NCFWF^T$. It is then possible to correct the FWF water path of each pixel.

<figure><img src="image_5.png" /><figcaption>**Figure 2.** FWF analysis of site #1 (see location in Figure 1f): (a) Initial FWF relative intensity $DN$ versus range $R$ presented between $R_b$ and $R_e$ bounds with a narrow filter (1) and large 3-m FWHM Gaussian low pass (2) allowing the evaluation of the offset between sea bottom and last echo without noise; (b) unscaled shape of the 532-nm emission pulse (1) compared to $dNCFWF^T$ (2) prior water index correction and $d(NCFWF^T)^{-1}$ (3) after water index correction; (c) wide filtered dFWF (1) and threshold level of relative intensity (2) used to delimitate the bounds $R'_b$ and $R'_e$ of CFWF calculation, dddNCFWF (3) used for echo determination in the range interval $R'_b-R'_e$; (d) transposition of the unscaled shape of 532 nm emission pulse (1) shown at arbitrary range, last 1064 nm discrete echo of water surface elevation (2), last dddNCFWF echo range with water index correction (3), transposition of $NCFWF^T$ without 1 (4) and with 1.34 (5) water index correction and its bottom range (6) at an offset of 2.8 m falling onto its corresponding last dddNCFWF echo.</figcaption></figure>
Following the Snell–Descartes law, as in [13] for example, light speed slows down from air $c_a$ to water $c_w$ and the incident angle of a laser beam in air $a_a$ rotates to $a_w$, according to their respective refractive indices $n_a$ and $n_w$. Since index $n_a$ is equal to 1 the 532 nm laser beam in water slows down proportionally to $1/n_w$ and the incident angle in water becomes $arcsin(sin a_a/n_w)$. Using a water refractive index of 1.34 (for a typical 35 psu salinity at 10 °C) for sea water in the La Baule area, a light round trip recorded in 1 ns in the air, defining its range resolution at 0.150 m, slows down 0.746 times to a new height resolution of 0.112 m in water. For example, a maximum apparent water bottom depth of the study area rises from 7 m to 5.22 m after water index correction in Figure 3a. Considering a maximum incident angle in the air of 14° rotation to 10.4° in water in Figure 3a, the lateral shift of a water FWF bottom would be smaller than 0.32 m, a third of a pixel width, and the vertical shift would be 0.07 m corresponding to half the FWF range resolution in the air. Both corrections are therefore disregarded. However, in strip borders imaged with an incident angle of 14°, an undulation of water surface of ±10° would project the FWF with a vertical error of up to 0.108 m, near the FWF range resolution in water, whereas the maximum horizontal error of 0.57 m remains within a pixel resolution. A vertical correction is therefore necessary along strip borders. The ratio between the vertical range of the incident angle in air ($\alpha_i$) and its range in water ($R_w$), xa and xw respectively in Figure 3, is equal to the ratio between the cosine of incident angle in water $\alpha_w$ and the cosine of incident angle in air $\alpha_a$. Therefore, at a 1-m-wide pixel location with an incident angle $\alpha_a$ in air with a 28° FOV, the full vertical range water index correction becomes:

$$NCFWF^T(Rw_{536nm}) = Ra_{1064nm} + \frac{(NCFWF^T(Ra_{536nm}) - Ra_{1064nm})}{n_w} \times \frac{\cos(\alpha_w)}{\cos(\alpha_a)} \quad (2)$$

Much higher incident angles would require wave slope corrections as shown by [35], which also justifies our choice of 28° FOV.

<figure><img src="image_6.png" /><figcaption>**Figure 3.** Cube image focalized on a plane at elevation 0 m visualized on (x, z) section with a record of the incident angle $\alpha$: (a) Example of light path rotation with sea waves undulating from -10° to 10° on a maximum incident angle of laser pulse at 14° showing that all apparent ranges of -7 m fall into a single 1 m (x, y) pixel and 0.15 m (z) elevation in air and at -5 m in water with a 1.34 water index correction; thick line caption: apparent range (R) in air is orange, shortened R in water is green and light blue with refraction, water surface is blue; thin light blue lines are refracted LiDAR ranges due to -10° to 10° water surface undulations; (b) (x, y) plane of the scanline geometry giving a row spacing along track of 0.5 m at nadir to 1 m on side strip image with a minimum point density of 1 m⁻¹ along track. Because of the 3 × 3 pixel median used for noise reduction, shows a three-pixel segment (a) and three rows in map (b), the effective resolution is 3 m × 3 m in both row and line directions.</figcaption></figure>
A full correction of the refraction angle would require a rotation of the raytracing vector at $Ra_{1064nm}$ and therefore a new projection of the FWF. But the resulting translation remaining within a pixel area would not improve the final results. This is particularly the case along strip borders where the point density goes down to $1 \text{ m}^{-1}$ as shown in Figure 3b. Moreover the $3 \times 3$ pixel median used to reduce the noise retaining only 1 best out of 9 neighboring pixels the effective resolution is $9 \text{ m}^2$ and a translation down to 6 m depth is not detectable. Each pixel of the image cube contains only one tight scan angle preserving the signal of interaction of the laser beam with a material. Such data arrangement is not compatible with light path line braking at varying water levels. This would require returning to the projection step to again process the same data with all the necessary parameters (pixel air/water elevation, slope and incident angle rotation) for a full refraction correction and data voxelization, for tomography purposes, which is out of the scope of this work.

The $NCFWF^T$ curves after water correction, can be transposed back into the original orthonormal plot in $(NCFWF^T)^{-1}$ at the single range resolution in air regardless of the water surface elevation. In practice a multiplication by 10,000 of both normalized $NCFWF$ and $(NCFWF^T)^{-1}$ functions is used to facilitate the comparison between their derivatives and original FWF at comparable 16 bits scale also allowing final image storage as an integer although all calculus are done in floats. This step is also used to cutoff last digits which ensures closure of the $NCFWF$ to 10,000 avoiding annoying noise. A gradient $\{-1, 0, +1\}$, followed by two low-passes $\{0.11, 0.22, 0.34, 0.22, 0.11\}$ for $NCFWF$ and $(NCFWF^T)^{-1}$ give the first smoothed derivative of both functions respectively $dNCFWF$ (without) and $d(NCFWF^T)^{-1}$ (with water index correction) in Figure 2b. Because of the water index correction $d(NCFWF^T)^{-1}$ displays a baseline tilted-up on the deep-water side. A simple shift proportional to the difference of relative intensity between beginning and ending of this function tilts back the baseline on its original state prior to its display in Figure 2b. As the water correction shrank the extent of $d(NCFWF^T)^{-1}$ and reduced the separability between echoes (Figure 2b), we used $dNCFWF$ for the echo calculations. The second derivative of $NCFWF$ given by the gradient $\{+1, 0, -1\}$, also followed by two low-passes $\{0.11, 0.22, 0.34, 0.22, 0.11\}$, application to $dNCFWF$ gives $ddNCFWF$. The third derivative of $NCFWF$ with the gradient $\{-1, 0, +1\}$ and two $\{0.11, 0.22, 0.34, 0.22, 0.11\}$ low-passes application to $ddNCFWF$ gives $dddNCFWF$ whose local maximum peaks, greater than a threshold chosen by the operator, provides a selection of echoes. Although rarely used the same procedure applied to $d(NCFWF^T)^{-1}$ gives $ddd(NCFWF^T)^{-1}$. The gradient/smoothing operations workflow is summarized in Figure 4.

<figure><img src="image_7.png" /><figcaption>**Figure 4.** Gradient/smoothing operations workflow. The bathymetry can be reached by direct calculation of echoes with $ddFWF$ on the initial raytracing image, by the first method with echoes derived from $dddNCFWF$, by the second method without echo based on bottom of $NCFWF^T$ or, in the case of clear water with well-formed echoes, by the third method with echoes derived from $(NCFWF^T)^{-1}$.</figcaption></figure>
A gradient {-1, 0, +1}, followed by 2 {0.11, 0.22, 0.34, 0.22, 0.11} low-passes applied to NCCFWF can also give the dNCCFWF useful for FWF damping function analysis as shown by [22].

Site #2 (Figure 5a) illustrates the effect of local FWF signal drop (such as electronic ringing effect) on the detection of small echoes. Site #3 (Figure 5b) is a typical case of FWF without local maximum at its bottom allowing the detection of echo. It is therefore important in such cases to determine the bathymetry by looking at the range of the final NCFWF<sup>T</sup> damping function in Figures 2d and 5b. It displays an offset to the effective range of bottom level echo proportional to the length of the laser pulse emission. The difference between the top and bottom of this FWF emission example is equal to 2 m with a maximum much closer to its top because of its asymmetry as shown in Figure 2b. But all echoes returned display a larger function because of their interaction with different medium on the ground, in the air and water. The offset between the bottom of NCFWF<sup>T</sup> and the last echo must be determined for each study. It is 2.8 m for a sandy bottom in both sites #1 and #3, Figure 2a,d and Figure 5b. Its full determination is part of Section 3 using GPS ground control points.

<figure><img src="image_8.png" /><figcaption>**Figure 5.** FWF analysis in two special cases #2 and #3 (see location in Figure 1f): (a) Narrow (1) and wide (2) low-pass of the FWF in site #2; (b) Narrow (1) and wide (2) low-pass of the FWF in site #3.</figcaption></figure>

The depth of the bathymetry is therefore available either by echo determination on dddNCFWF when a local maximum, eventually enhanced by successive applications of short high pass and long low pass filters, can be detected or at the base of NCFWF<sup>T</sup> when no local maximum is available. Consequently, two bathymetric methods are available: (1) the last echo can be deduced from the range of the last local maximum of dddNCFWF followed by a water index correction; (2) the bottom of the water index corrected NCFWF<sup>T</sup> can be used to infer the range of the last unshaped echo with the addition of an offset determined by ground control points.

### 2.2.3. GPS Ground Acquisition

A GPS acquisition campaign was carried out at low tide in the western part of the bay in order to obtain the foreshore topography (Figure 6). More than 90 control points were levelled between the 0 Chart Datum closely corresponding to the lowest astronomical tide (see Table 1) to the top of the dry beach, passing over foreshore forms such as sand banks.

<figure><img src="image_9.png" /><figcaption>**Figure 6.** Field sampling of GPS heights by Martin Juigner: example of site #40 showing emerging rocks covered with brown algae.</figcaption></figure>
The reference altimetry for the metric units of this campaign is NGF-IGN69. In order to establish seamless data between topographic and bathymetric references, we used the calculations of the differences between the bathymetric reference system and the topographic reference system listed in Table 1 for two calculation points surrounding the bay: Pouliguen to the west and Pornichet to the east.

<table><caption>**Table 1.** Altimetric references in meters [36].</caption><thead><tr><td>Meter (Chart Datum and NGF-IGN69)</td><td>Le Pouliguen</td><td>Pornichet</td></tr></thead><tbody><tr><td>Highest Astronomical Tide (CD ref)</td><td>6.15</td><td>6.22</td></tr><tr><td>Mean Sea Level (CD ref)</td><td>3.31</td><td>3.41</td></tr><tr><td>Lowest Astronomical Tide (CD ref)</td><td>0.04</td><td>-0.01</td></tr><tr><td>Chart Datum in NGF-IGN69 ref</td><td>-2.917</td><td>-2.927</td></tr></tbody></table>

### 2.2.4. Multibeam Data Acquisition

Bathymetric data were acquired in May 2018, during the Hoopla18 geophysical survey with R/V Haliotis (Ifremer/Génavir). It is equipped with a hull mounted interferometric sonar (GeoAcoustics, GeoSwath, 250 kHz) and a chirp seismic profiler (1.8 to 5.3 kHz), allowing simultaneous acquisition of sonar bathymetry, sonar imagery and seismic profiles. Information from an Aquarius Thales GPS and an RTK beacon on the coast, served to position the data with centimeter precision and appropriate tidal corrections. All the bathymetric data were processed by [37] with Caraïbes software [38] from [39]. Two profiles were acquired from 5 m to 4.5 m of water depths, in the east of the Bay, perpendicular to the beach, in front of the harbor. These two profiles gave calibrated measurements both in bathymetry and position. Thus, they may be used as calibration data for further acquisition. The seismic data confirm that the sea bottom is covered by sandy formations.

## 3. Results

### 3.1. La Baule Bay

The ray tracing FWF image of La Baule bay mosaic (Figure 7a) has 20 strips made of 53 segments. A quick view of the ray tracing FWF image can be made with three channels at the elevation -1 m, -4 m and -7 m NGF respectively in red, green and blue of the composite image in Figure 7a. The mean, max and min color composite image, Figure 7b, similar to Figure 1c,d, displays different false colors. This difference is from numerous baseline drops modifying the distribution of the minimum values, highlighted in blue, all pixels free of intensity drops.

The thickness of the FWF range characterizes the laser beam diffusion or scattering in a pixel. It varies from the bottom of the water to its surface at sea and from the ground to the top of trees on land (Figure 7c). When transformed into NCFWF this diffusion range includes the laser pulse thickness which explains why all pixels display a minimum offset corresponding to the laser pulse length of interaction with a flat surface. In Figure 7c this offset after NCFWF transformation is greater than 5 m, with rare exceptions, justifying a detection limit to 4.5 m. In comparison with the 2-m effective length of laser pulse emission example in Figure 2b it is clear that such minimum length of FWF on beach and foreshore comes from the interaction of that laser beam with dry and wet sand respectively. In order to check any effect of incident angle this data is presented in Figure 7d as a control.

Large FWF signal drops along beaches, as exemplified in Figure 5a, and the weakness of echoes in turbidity area limit the bathymetry mapping by ddFWF in Figure 8a. The echo detection is largely extended when working with dddNCFWF in Figure 8b but remains impossible when dddNCFWF display no local maximum. Only bottom of NCFWF<sup>T</sup> extends the bathymetric mapping throughput the bay in Figure 8c with a noticeable exception in the south where turbidity and wind wrinkles hide some of it. Considering the offset thickness observed in Figure 2a a shift of 2.8 m was applied to the NCFWF<sup>T</sup> bottom in Figure 8c.
<figure><img src="image_10.png" /><figcaption>**Figure 7.** FWF: (a) color composite image of elevation -1 m, -4 m and -7 m in red, green and blue respective channels, (b) minimum, maximum and mean intensities color composite image of 532 nm FWF, (c) row thickness of NCFWF ranges, (d) 532 nm laser beam incident angle.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>**Figure 8.** FWF extraction of echoes: (a) from ddFWF, (b) from dddNCFWF; (c) bottom of the NCFWF<sup>T</sup> shifted by 2.80 m and (d) same data with inconsistent data masking Figure 9. Echoes and NCFWF<sup>T</sup> bottom elevations are all shown with full range water index correction.</figcaption></figure>

A water column greater than 4 m being inconsistent with a bathymetry lower than 4 m, both criteria are combined with pixels having last 532 nm echoes between 0 m and 5 m to build the mask in Figure 9a. One can also use the dNCCFWF for a more complete analysis of the laser beam light path in Figure 9b. It is a color composite image showing large downward diffusion, down to 5 m depth in sea water (in blue) but also upward diffusion, up to +1 m, where the bottom displays a stronger echo than the water surface (in orange), whereas homogeneous less diffusive dry sand appears dark brown at the top of beaches in Figure 9b.
<figure><img src="image_12.png" /><figcaption>**Figure 9.** Complementary data: (a) inconsistent data masking (black area) by a combination of 3 masks: in red for bottom NCFWF$^T$ > −4 m, green for last green echoes > 0 m and < 5 m, blue for FWF water height > 4 m; (b) dNCCFWF color composite image of negative downward FWF diffusion in blue and green and positive upward diffusion in red from the FWF maximum (or main echo).</figcaption></figure>

## 3.2. Low Tide GPS Ground Truthing

The 92 GPS ground control sites, located on the Benoit beach area orthophoto Figure 10a, were sampled (see Figure 6) for comparison with the LiDAR range results in Figure 11.

<figure><img src="image_13.png" /><figcaption>**Figure 10.** Benoit beach GPS ground sampling: (a) orthophoto showing finger shapes of turbidity current between water layers w1 to w3, (b) FWF color composite image with −6 m, −3.2 m, −1 m in red, green, and blue channels respectively, (c) last echoes of dddNCFWF, A-B is the location of the range slices presented in Figure 12. (d) intensities of the last dddNCFWF echoes and (e) bottom of NCFWF$^T$ with a 2.8 m offset; pixels 1 and 2 are presented in Figure 13. The stair step range in (c) is due to the FWF range resolution of 0.15 m. Both dddNCFWF and NCFWF$^T$ maps are elevations shown in (c,e) with full range water index correction.</figcaption></figure>
GPS-derived z elevations, (1) in Figure 11a, draw a diagonal with a slope of 1 in the *R* versus *z* plot in Figure 11a like in the *z* versus *R* plot of Figure 11b while Figure 11c displays the spatial distribution of *R* and *z* along the pathway of the 92 ground control sites. The LMS 532 nm discrete echoes, (2) in Figure 11, provided by LMS are well aligned on a regression line with a correlation slope equal to 1.08 and reciprocally with a slope of 0.91 in Figure 11b. The 1064 nm last discrete echoes, (3) in Figure 11a provide the sea level (Figure 11a,c) at different tides during Day 1 and Day 2. Despite this tide level variation from Day 1 to Day 2 all first method water index corrected echoes, (4) in Figure 11, deduced from dddNCFWF are aligned on a straight line in both Figure 11a,b giving a slope correlation of 0.96 with *z* in b which will be used to calibrate all subsequent first method results. For comparison with the second method, the bottom of water corrected NCFWF<sup>T</sup>, (5) in Figure 11, is also aligned with a slope of 0.95 characterized by an offset of 2.8 m in Figure 11b, which will be used for calibration of second method results. It also confirms the observation made in Figure 2a. The addition of this constant 2.8 m offset to the NCFWF<sup>T</sup> bottom alone, (6) in Figure 11, allows the retrieval of elevation *z*. The plot along site numbers in Figure 11c draws a type of bathymetric profile showing that LMS deduced echoes (2) and dddNCFWF echoes (4) given by the first method are better aligned on GPS control point *z* than echoes of the second method (6) derived from NCFWF<sup>T</sup> bottom forming a less smooth profile.

<figure><img src="image_14.png" /><figcaption>**Figure 11.** Comparison between *z* elevation determined by ground GPS measurements and various LiDAR ranges, with water index correction when applicable, obtained along 92 sites: (a) plot of *R* results versus *z*; (b) plot of *z* versus *R* giving the correction coefficients; (c) display of all elevations (*z*) along the 92 site numbers with 2 and 9 weighted by 0.96. Common plot caption: (1) GPS derived *z* elevation, (2) discrete 532 nm last echoes provided by LMS, (3) discrete 1064 nm first echoes, (4) bottom of NCFWF<sup>T</sup> in water, (5) bottom of water NCFWF<sup>T</sup> in water plus 2.8 m offset, and (6) dddNCFWF last echoes with water index correction. The GPS-derived elevations *z* are plotted on *R* axis with 0° incident angle.</figcaption></figure>

Turbidity at the Benoit beach area is sufficiently low to see bottom rocks which are partially covered in brown algae (Figure 6) in the south between sites #30 and #40 (Figure 10). The apparent reflectance appears greater than 20% because of low turbid water column diffusion. Due to tide
currents at water bottom, turbidity currents appear as fingers near site #60 in Figure 10. At least three
types of water color, from w1 to w3 in Figure 10a, form digitation in each other at their interface.
This indicates the presence of water layering which is also well observed in the color composite image
of the corresponding FWF in Figure 10b. This variation of turbidity has no impact on the range of
dddNCFWF echoes in Figure 10c and only changes their intensity in Figure 10d. Turbidity currents
appear to be associated to the accretion slope of two sand banks slowly moving towards the north and
forming a depression between them. It is clearly visible in Figure 10c but partially masked by currents
of turbidity in the NCFWF<sup>T</sup> bottom with 2.8-m offset of the second method in Figure 10e. This latter
method seems to be more sensitive to the highest turbidity in contact with sand remobilization layer at
the sea bottom whereas the dddNCFWF echoes pick-up the whole shape of the interbank depression.

The range slice of FWF allows better understanding of why some echoes are missing. In Figure 12a
there is a continuous dropping of signal without local maximum preventing the detection of echoes.
Figure 12b displays the dNCFWF range slice, whose bottom is equal to the NCFWF<sup>T</sup> bottom prior
water index correction, showing that the stair step bottom detection of the dNCFWF function allows
the localization of the sea bottom with constant offset. Because of its continuous decreasing shape,
the dddNCFWF cannot display local maximum along the range slice in Figure 12c and echo detection
fails locally (see arrow). The ddd(NCFWF<sup>T</sup>)<sup>-1</sup> function after water index correction gives the right
bathymetry in Figure 12d but its combination with water level detection by 1064 nm discrete echoes
reinforces the noise and more sea bottom is missing (compare d with c) which explains why we did no
choose to process ddd(NCFWF<sup>T</sup>)<sup>-1</sup> data in the following work.

<figure><img src="image_15.png" /><figcaption>**Figure 12.** Range slices of: (a) FWF, (b) dNCFWF, (c) dddNCFWF and (d) full range water index corrected ddd(NCFWF^T)^{-1}, all in relative intensities (DN) of Profiles A-B from Figure 10. The arrow in (a) points towards continuous FWF without bottom echo in (b–d).</figcaption></figure>

The turbid water layering observed in Figure 10 produces interfaces between them. As discussed
in the introduction, both FWF echoes at sea surface and sea bottom are linked to each other by
intermediate values due to the intensity of laser backscattering along the column of water as shown
in the two examples in Figure 13. However, the example presented in Figure 13 was obtained with a softer processing using only one low-pass allowing the visualization of intermediate layer while it is not visible in Figure 12. Its effective detection would also need the detection of local maximum in the negative part of the dddNCFWF function. Nevertheless, this proves that turbid water can be made of a succession of echoes progressively masking the last bottom echo in a single continuum by filling the gap between the first and last echo.

<figure><img src="image_16.png" /><figcaption>**Figure 13.** Profiles of FWF and dddNCFWF relative intensities (DN) of two pixels sampled in Figure 10 with only 1 low-pass application allowing the visualization of water layering instead of two low-passes used to minimize them.</figcaption></figure>

So, without a priori knowledge of the composition of the sea bottom and the sediment mobilized in the water column, we were able to detect the right bathymetry of the sea bottom as shown in Figure 11 with a calibration to GPS ground control points acquired at low tide (Figure 6) using: (1) 0.96 gain correction on dddNCFWF echoes with the first method and (2) 2.8 m offset on NCFWF<sup>T</sup> bottom with the second method of bathymetry mapping.

### 3.3. Multibeam Bathymetry Truthing

The multibeam data were reprocessed at FWF spatial resolution (1 m × 1 m) for a more rigorous comparison with dddNCFWF echoes. Both lines of acquisition are mosaicked in Figure 14a with color coding on a −5.5 m to −5 m range highlighting the presence of a narrow sandbank set (see arrow) while Figure 14b displays the same data in a range of −6 m to 6 m facilitating the comparison with Figure 8. The first method based on dddNCFWF echo detection fails to map the bathymetry in the eastern part of Figure 14c because of higher turbidity and superficial wind effects.

On the contrary, Figure 14d, the second method based on 2.8 m offset of NCFWF<sup>T</sup> bottom can map all the bathymetry in the full area despite a stronger noise surprisingly reinforcing the contrast of sandbanks on the sea bed.

A profile A-B drawn on each data from Figure 14 allows the quantification of the deviation between calculated bathymetry with LiDAR data and the multibeam echo sounder bathymetry in Figure 15. The multibeam profile, (1) Figure 15, confirms the general flat distribution of the bathymetry in the eastern part of La Baule bay. Most of it is around −5 m NGF with a steeper slope next the beach and Pornichet harbor (see location in Figure 1). The discrete 1064 nm last echoes, (2) Figure 15, give the sea level. Day 1 data displays a gentle slope from strip 8 to strip 10, marking the tide progression with time. Day 2 from strip 11 to strip 15 displays the same slope plus well-formed waves of oceanic swell oblique on the profile with a period length of ~100 m giving slopes at ±6°. The first method last echoes based on dddNCFWF, (3) Figure 15, provides a smooth bathymetry but unfortunately with blind areas where the water surface seems to stop the laser beam propagation (see local profile jump to the water surface level). The second method derived from bottom of the water index corrected NCFWF<sup>T</sup>, (4) Figure 15, is ~2.8 m below the reference level (1) and the application of the constant offset of 2.8 m, (5) Figure 15, provides a good approximation of the bathymetry despite a stringer level of
noise with a 0.33-m standard deviation to multibeam data. Such bathymetry is therefore provided by
the second method within a 2-sigma interval of 0.66 m around the effective bottom.

<figure><img src="image_17.png" /><figcaption>**Figure 14.** Multibeam comparison: (a) Multibeam data showing sandbank orientation (arrow), (b) same multibeam data at scale of FWF data, (c) first method dddNCFWF last echoes with the 0.96 correction using Benoit beach GPS control points, (d) second method bottom of NCFWF<sup>T</sup> range with 2.8 m offset.</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>**Figure 15.** Comparison between bathymetric profile A-B across strip 8 to 15 from Figure 14 of: (a) (1) multibeam bathymetry, (2) 1064 nm discrete last echoes, (3) last dddNCFWF echo time 0.96, (4) bottom of water index corrected NCFWF<sup>T</sup>, (5) same NCFWF<sup>T</sup> bottom with 2.8 m offset, (6) laser beam incident angle; (b) FWF range slice and (c) dddNCFWF range slice (without water index correction) along the same profile A-B.</figcaption></figure>
The FWF range slice along profile A-B of Figure 14 presented in Figure 15b clearly shows an increasing intensity near the bottom. The long recording of FWF down to -20 m (before water index correction) was important for the detection of the right baseline of each pixel. The sea bottom is well defined on the dddNCFWF range slice in Figure 15c but the missing part still justifies the use of the second method.

## 4. Discussion

The method presented in this work uses a straight ray tracing procedure to build a mosaic of FWF on a 1-m-wide pixel grid with a vertical range resolution of 0.15 m along a laser pulse light path with preservation of incident angle information in each pixel. Because of the narrow scanning angle used, the light path in water was successfully corrected along a 0.15 m range direction in the initial FWF data cube geometry without translation between pixels. This step could have been improved by considering the exact orientation of the water interface with air. However, the strong hot spot effects observed in water and most wet areas, such as sediment ripple marks exposed at low tide, suggest that capillary waves [40] can play an important role as shown also by [21]. A full characterization would therefore require a lot more work and a much higher resolution note compatible with the monitoring of hundreds of shoreline kilometers. At a scale of 1-m-wide pixel for a bathymetry down to 5-m depth NGF, the fast procedure used seems appropriate. A 3-m wide pixel for 10 m depth bathymetry is also possible for large shoreline monitoring. Above that level marine applications are more time efficient. La Baule bay data was acquired on two different days because of temporary technical problems in 2017. A new campaign in 2018 covering 16.7 km² was done in a single day in less than 3 h.

As for [22] all figures presented in this work use nearest pixel sampling without pixel aggregation and all calculations were done at the 1 m pixel resolution to prevent any FWF stretching on the slope. The pixels in Figure 1 are 1 m square, but only 1 out of 100 are shown. Aggregation of pixels are only possible on final maps after range and echo extraction. The method relying on the FWF shape analysis, any distortion of it can move the range of echoes.

Radiative transfer procedures require a lot of prior knowledge of the sediment reflectance and attenuation coefficient of the water loaded with various particles and molecules as is commonly done with hyperspectral studies [14]. The proposed method removes most of these parameters through the normalization of the FWF cumulate whose derivative is only proportional to variation of interface elevation returning a signal whatever the composition of the bottom and the water column. It only considers the salinity of the water for the determination of the best water index allowing to slow down the light speed round-trip and determines the right level of sea bottom. Dark and bright targets look identical but of course noise level becomes an issue on darker medium. This is particularly the case when no echo appears at sea bottom because of multiple tiny returns due to the accumulation of intermediate turbid water layer returns as seen in Figures 12 and 13. We limit this effect by using a large Gaussian convolution filter to determine the bounds $R'_b$ and $R'_e$ of a useful segment of FWF (Figure 2c). A series of smoothing kernels are then used alternately with high-pass calculation of derivatives summarized in Figure 4.

Figure 16 summarizes the FWF typology highlighted in this work. In Figure 16a, all methods are applicable in clear water with sharp double-edged echoes at sea bottom. In Figure 16b, the first method based on dddNCFWF allows the detection of weaker double-edged echoes but the second method based of NCFWF$^T$ must be used for the detection of single edged echoes at sea bottom, Figure 16c. The blind areas (Figure 16d,e) correspond to very long damping function where no edge marks the end of FWF signal. The inconsistency between FWF potential water thickness and the bathymetry reached by the bottom of NCFWF$^T$ was used to delineate such blind areas in Figure 9a and combined with binary its dilation to mask the full blind areas in Figure 17.

From the median applied on 3 × 3 pixel boxes, to remove bad returns in air, all the processes proposed in the present method work at a pixel scale without any consideration for the neighboring values and without any interpolation. This is particularly important to avoid laser pulse stretching
on the slope. It is however possible to look at a large neighborhood in the final map of sea bottom to produce smoothed bathymetry. The following map, Figure 17, is an aggregation of 16 pixels smoothing and reducing its spatial resolution to 4 m wide pixels. A 0.5 m noise appears along the strip border on data from Day 2 whereas Day 1 data are free of noise in Figure 17. The first day is characterized by a low tide and flat sea surface (Figures 11 and 16) whereas the second day displays stronger waves and a high tide as shown in the profile in Figure 15. This strip border noise mainly occurs for incident angles greater than 14° (Figure 7d) during high tide with well-formed sea waves. At low tide the conditions of observation were a lot greater with less noise and lower sea waves. A more careful adjustment of the airborne mission from day to day with tide conditions would save a lot of time in data processing. The better the data acquisition the faster and better the result. Early morning without thermal wind conditions is also recommended.

<figure><img src="image_19.png" /><figcaption>**Figure 16.** FWF typology: (a) clear water with distinct double-edged (2 arrows) echoes at bottom; (b) light turbid water with double-edged echoes connected by intermediate echoes; (c) moderate turbid water with continuous top to bottom merged intermediate echoes with a single edge (1 arrow) at the bottom; (d) turbid water with beveled end of FWF without edge (x) at bottom; (e) high turbidity with narrow echoes at sea surface.</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>**Figure 17.** NCFWF<sup>T</sup> bottom bathymetry map with aggregation of 16 pixels after calculation of a three-pixel size median on the final map; HF are hidden faults driving the morphology in Le Pouliguen—La Baule Benoit beach area.</figcaption></figure>
Despite a stronger noise level than echoes obtained with the dddNCFWF first method, the NCFWF<sup>T</sup> bathymetry of the second method provides a wider coverage and gives the best morphology map of sea bottom for sand bank tracking in the low turbidity bay of La Baule (Figure 17). In Le Pouliguen area and Plage Benoit, hidden faults (HF in Figure 17) delimit underlying tilted stair step granite blocks. Wide, high sandbanks are growing along both hidden faults forming weak crests along the shoreline in Le Pouliguen-Plage Benoit while low narrow sandbanks perpendicular to the shoreline characterize the main zone of erosion of La Baule-Pornichet area. Wide sandbanks are also characterized by asymmetric profiles as shown by the dddNCFWF range slice in Figure 18a and its corresponding dddNCFWF last echo profile in Figure 18b with a long slope of accretion (arrows) gently dipping southward and a short steep slope dipping northward (see also Figure 10). Most remarkable rocks and sea bed morphology are also clearly identifiable in Figure 17 map as in the dddNCFWF range slice along profile C-D in Figure 18c. The blind areas correspond to a 6 m depth artifact along the C-D profile.

<figure><img src="image_21.png" /><figcaption>**Figure 18.** NCFWF<sup>T</sup> bottom bathymetry: (a) dddNCFWF range slice in air and (b) its corresponding last echoes with water index correction of the A-B profile (Figure 17) showing the wide sandbank asymmetry; (c) transversal profile C-D (Figure 17) of dddNCFWF range slice in air showing intermediate layers in the blind areas.</figcaption></figure>

**5. Conclusions**

The definition of a water column as an image cube of 1-m-wide pixels with 0.15-m range resolution spectra of a 532-nm laser beam relative intensity returned by multiple turbidity layers allowed mapping the bathymetry of sea bottom with FWF records while 1064 nm discrete echo records ensure the determination of water surface. A maximum scan angle of 14° (28° FOV) project all FWF in one pixel when the water column is lower than 10 m and the wave slope is lower than 10°. This in turn reduces the water index correction to a height resolution of 0.112 m below water surface varying with the wave elevation of each pixel. The increasing turbidity progressively adds new returns forming a single continuum from top to bottom preventing the formation of a complete echo at the bottom. A wide Gaussian filter helped to delineate the range limits of the meaningful FWF from the blank noise of a baseline which must be determined for each pixel.

Because of the wide variation of intensities ranging from strong specular reflections on superficial waves down to weak reflections of bottom rocks covered by brown algae, a drastic simplification is done by the normalization of a cumulative FWF giving the same weight to bright and dark FWF pixels. Each pixel is then a spectrum in range of relative intensities proportional to the distribution of returns along the water column at sea. The smallest FWF range is given by one pixel of sand diffusing 532 nm and 1064 nm laser-beams in all directions and returning a narrow echo for each laser wavelength. The typical echo formed on wet sand in shallow clear water displays a downward damping function of 2.8 m whereas it is only 2 m for the emission pulse which confirms that its diffusion at water/sand interface tends to enlarge the FWF. At the water surface the hotspot specular effects mainly increase the intensity without enlarging the FWF too much which can however locally block the observation of sea bottom. In water, the multiplication of laser-beam backscattering on intermediate interfaces of turbid water layers tends to connect top and bottom echoes in a single damping function. Thus, an FWF ending should display a final 2.8-m-wide stair step edge at the end of its damping function on
a real bottom. The edge detection at a distance over the bottom can also occur on an unexpected noise layer found at 6 m below sea level and masking the bottom in Figure 18c.

The 2.8 m offset between last echo and the bottom of the FWF damping function given by NCFWF<sup>T</sup> was confirmed in the Benoit beach area by a regression line of correlation with GPS derived elevation. This offset value was additionally validated with multibeam echo sounder data down to 5.5 m depth NGF. The proposed method of bathymetry is therefore semi-empirical since it necessitates GPS ground control points at low tide and independent multibeam echo sounder to confirm the offset between NCFWF<sup>T</sup> bottom and effective sea bed elevation.

The precision of the first method dddNCFWF echo's elevation is given by the standard deviation with multibeam elevations of 0.13 m in the Pornichet area, below the 0.15-m range resolution of FWF recording. This standard deviation increases to 0.33 m in the same area with the second method NCFWF<sup>T</sup> bottom which means that results are within a layer of 0.66 m.

From the pixel aggregation in Figure 17 and the initial median filter for atmospheric noise removal, all calculi were done exclusively along the incident angle vector in the FWF range dimension without considering any neighboring data. Of course, classical data cloud interpolation between neighbors could have helped to improve the bathymetry map. Other tomographic methods based on voxel projections or cube image processing could also be used for the detection of the lateral extension of multiple interfaces between various water layers. This will be done in a further study combing synchronous hyperspectral images digitized in the same incident angle directions.

The main goal of this study was to show that a pure spectral analysis of FWF damping function, exclusively along tight incident angle and range dimension, can map a bathymetry with (first method) or without (second method) echo detection in moderately turbid bays. All normalizations are reversible at any time in the proposed workflow and diagnosis of bottom are also possible. The interest of the method is evidenced by the application to a case study at La Baule with the monitoring of wide along beach sandbanks presenting turbid current on accretion slopes in the Plage Benoit area. It is also attested by the concentration of narrow and low radial sandbanks in the erosive Pornichet area along the direction of tidal currents. The detection of turbid water layering needs to be further investigated with synchronous sampling in the water column and periodic artifact removal.

**Supplementary Materials:** The following are available online at http://www.mdpi.com/2072-4292/11/2/117/s1, S1_High_resolve_figs: high resolution Figures.

**Author Contributions:** P.L. and M.G. conceived the study and analyzed the data; M.R and A.B. contributed to the validation; P.L. wrote the paper.

**Funding:** This research was supported by the Region Pays de la Loire with funding of the RS2E-OSUNA and OR2C programs. The Nantes Rennes LiDAR platform was funded by the Region Pays de la Loire and the Region Bretagne with the European Regional Development Fund (ERDF).

**Acknowledgments:** We would like to thank William Gentile, Cyril Michon and Emmanuel Gouraud from GEOFIT Company. Jean-Jérôme Houdaille from PIXAIR Company, Anca Dobrinescu, Paul La Rocque, Mark Topping and Michel Stanier, from Teledyne Optech Company. Martin Juigner from the LETG Nantes for GPS ground control points and Martin Sanchez from the LPG Nantes for turbidity data in La Baule. Particular thanks to Dimitri Lague from Geosciences Rennes and Erwan Le Menn from the LPG Nantes, for their respective management and contribution to the Nantes Rennes LiDAR platform and Laurence Huber Moy from the LETG in Rennes for her help in creating of the shared platform between Nantes and Rennes. We would also like to thank the reviewers who helped us to better focus this publication on methodology and main results.

**Conflicts of Interest:** The authors declare no conflict of interest.

## Appendix A

Because of the FOV set to 28°, the swath or strip image width was 199.5 m on the ground with a nominal footprint of 0.49 m along-track and 0.29 m cross-track. Because of the tilt of the 532-nm laser beam and because of a sidelap set to 30% to prevent empty space between strips, the effective incident angle range was from 4 to 20°. With an aircraft speed of 201.7 km/h the down track footprint varied from 0.49 m at nadir to 0.98 m on sideways defining a minimum resolution size ensuring complete
ground coverage without pixel gaps. With a laser pulse frequency set between 50 and 100 kHz (below 175 kHz for FWF optimal recording), the overall mean point density was 3.52 to 7.05 points/m². The FWF was recorded with a maximum length of 60 m, at 1 GHz frequency, resulting in a vertical resolution of 0.15 m. The roll compensation was also active to minimize geometric distortion during flights. The Optech LiDAR Mapping Suite (LMS) [23] combines a Global Positioning System (GPS) and Inertial Measurement Unit (IMU) associated to the LiDAR to provide a georeferenced point cloud and associated FWF with strip optimization. To compute the complete point cloud, LMS uses an analogous approach to photogrammetric block adjustment. On overlap areas between each flight line, planes (typically roofs) are extracted and used in a least squares adjustment model provided by LMS [23]. The trajectory accuracy provided by GEOFIT Company was better than 0.15 m in planimetry and 0.08 m on elevation.

## References

1. Del Rio, L.; Gracia, F.J.; Benavente, J. Shoreline change patterns in sandy coasts. A case study in SW Spain. *Geomorphology* **2012**, *196*, 252–266. [CrossRef]

2. Jackson, C.W., Jr.; Alexander, C.R.; Bush, D.M. Application of the AMBUR R package for spatio-temporal analysis of shoreline change: Jekyll Island, Georgia, USA. *Comput. Geosci.* **2012**, *41*, 199–207. [CrossRef]

3. Moussaid, J.; Ait Fora, A.; Zourarah, B.; Maanan, M.; Maanan, M. Using automatic computation to analyze the rate of shoreline change on the Kenitra coast, Morocco. *Ocean. Eng.* **2015**, *102*, 71–77. [CrossRef]

4. Dada, O.A.; Li, G.; Qiao, L.; Ding, D.; Ma, Y.; Xu, J. Seasonal shoreline behaviors along the arcuate Niger Delta coast: Complex interaction between fluvial and marine processes. *Cont. Shelf Res.* **2016**, *122*, 51–67. [CrossRef]

5. Mann, T.; Westphal, H. Multi-decadal shoreline changes on Takú Atoll, Papua New Guinea: Observational evidence of early reef island recovery after the impact of storm waves. *Geomorphology* **2016**, *257*, 75–84. [CrossRef]

6. Boak, E.H.; Turner, I.L. Shoreline definition and detection: A review. *J. Coast. Res.* **2005**, *21*, 688e703. [CrossRef]

7. Juigner, M.; Robin, M.; Debaine, F.; Hélen, F. 2017 A generic index to assess the building exposure to shoreline retreat using box segmentation: Case study of the Pays de la Loire sandy coast (west of France). *Ocean Coast. Manag.* **2017**, *148*, 40–52. [CrossRef]

8. Le Mauff, B.; Juigner, M.; Ba, A.; Robin, M.; Launeau, P.; Fattal, P. Coastal monitoring solutions of the geomorphological response of beach-dune systems using multi-temporal LiDAR datasets (Vendée coast, France). *Geomorphology* **2018**, *304*, 121–140. [CrossRef]

9. Pastol, Y. Use of Airborne LIDAR Bathymetry for Coastal Hydrographic Surveying: The French Experience. *J. Coast. Res.* **2011**, *62*, 6–18. [CrossRef]

10. Guenther, G.C.; Brooks, M.W.; LaRocque, P.E. New Capabilities of the “SHOALS” Airborne Lidar Bathymeter. *Remote Sens. Environ.* **2000**, *73*, 247–255. [CrossRef]

11. Lillycrop, W.J.; Irish, J.L.; Parson, L.E. SHOALS system: Three years of operation with airborne lidar ba thymetry—Experiences, capability and technology advancements. *Sea Technol.* **1997**, *38*, 17–25.

12. Irish, J.L.; Lillycrop, W.J. Scanning laser mapping of the coastal zone: The SHOALS system. *ISPRS J. Photogramm. Remote Sens.* **1999**, *54*, 123–129. [CrossRef]

13. Wang, C.-K.; Philpot, W.D. Using airborne bathymetric lidar to detect bottom type variation in shallow waters. *Remote Sens. Environ.* **2007**, *106*, 123–135. [CrossRef]

14. Jay, S.; Guillaume, M.; Blanc-Talon, J. Underwater target detection with hyperspectral data: Solutions for both known and unknown water quality. *IEEE J. Sel. Top. Appl. Earth Obs. Remote Sens.* **2012**, *5*, 1213–1221. [CrossRef]

15. Hapke, B. Bidirectional reflectance spectroscopy: 1. Theory. *J. Geophys. Res.* **1981**, *86*, 3039–3054. [CrossRef]

16. Mallet, C.; Bretar, F. Full-waveform topographic lidar: State-of-art. *ISPRS J. Photogramm. Remote Sens.* **2009**, *64*, 1–16. [CrossRef]

17. Kashani, A.G.; Olsen, M.J.; Parrish, C.E.; Nicholas Wilson, N. A Review of LIDAR Radiometric Processing: From Ad Hoc Intensity Correction to Rigorous Radiometric Calibration. *Sensors* **2015**, *15*, 28099–28128. [CrossRef]
18. Wagner, W.; Ullrich, A.; Ducic, V.; Melzer, T.; Studnicka, N. Gaussian decomposition and calibration of a novel small-footprint full-waveform digitizing airborne laser scanner. *ISPRS J. Photogramm. Remote Sens.* **2006**, *60*, 100–112. [CrossRef]

19. Ullrich, A.; Pfennigbauer, M. Echo Digitization and Waveform Analysis in Airborne and Terrestrial Laser Scanning. In Proceedings of the Photogrammetric Week, Stuttgart, Germany, 5–9 September 2011; Fritsch, D., Ed.; Wichmann/VDE Verlag: Belin, Germany; Offenbach, Germany, 2011; Volume 11.

20. Höfle, B.; Pfeifer, N. Correction of laser scanning intensity data: Data and model-driven approaches. *ISPRS J. Photogramm. Remote Sens.* **2007**, *62*, 1415–1433. [CrossRef]

21. Westfeld, P.; Maas, H.-G.; Ritcher, K.; Weiß, R. Analysis and correction of ocean wave pattern induced systematic coordinate errors in airborne LiDAR bathymetry. *ISPRS J. Photogramm. Remote Sens.* **2017**, *128*, 314–325. [CrossRef]

22. Launeau, P.; Giraud, M.; Ba, A.; Moussaoui, S.; Robin, M.; Debaine, F.; Lague, D.; Le Menn, E. Full-Waveform LiDAR Pixel Analysis for Low-Growing Vegetation Mapping of Coastal Foredunes in Western France. *Remote Sens.* **2018**, *10*, 669. [CrossRef]

23. Lidar Mapping Suite—Professiona; Teledyne Optech: Vaughan, ON, Canada, 21 November 2018.

24. Persson, A.; Sderman, U.; Tpel, J.; Ahlberg, S. Visualization and analysis of full waveform airborne laser scanner data. In Proceedings of the ISPRS Workshop “Laser scanning 2005”, Enschede, The Netherlands, 12–14 September 2005; pp. 103–108.

25. Wang, H.; Glennie, C.; Prasad, S. Voxelization of full waveform LiDAR data for fusion with hyperspectral imagery. In Proceedings of the IEEE International Geoscience and Remote Sensing Symposium—IGARSS, Melbourne, Australia, 21–26 July 2013; pp. 3407–3410.

26. Wang, H.; Glennie, C. Fusion of waveform LiDAR data and hyperspectral imagery for land cover classification. *ISPRS J. Photogramm. Remote Sens.* **2015**, *108*, 1–11. [CrossRef]

27. Buddenbaum, H.; Seeling, S.; Hill, J. Fusion of full-waveform lidar and imaging spectroscopy remote sensing data for the characterization of forest stands. *Int. J. Remote Sens.* **2013**, *34*, 4511–4524. [CrossRef]

28. Hermosilla, T.; Ruiz, L.A.; Kazakova, A.N.; Coops, N.C.; Moskal, L.M. Estimation of forest structure and canopy fuel parameters from small-footprint full-waveform LiDAR data. *Int. J. Wildland Fire* **2014**, *23*, 224–233. [CrossRef]

29. Nie, S.; Wang, C.; Zeng, H.; Xi, X.; Li, G. Above-ground biomass estimation using airborne discrete-return and full-waveform LiDAR data in a coniferous forest. *Ecol. Indic.* **2017**, *78*, 221–228. [CrossRef]

30. LAS Specification Version 1.4—R13; American Society for Photogrammetry & Remote Sensing: Bethesda, MD, USA, 2013; 28p.

31. Jähne, B. *Digital Image Processing, Concepts, Algorithms and Scientific Applications*; Springer: Berlin, Germany, 1991; 383p.

32. Underwood, E.E. *Quantitative Stereology*; Addison-Wesley Publishing Company: Menlo Park, CA, USA; London, UK; Don Mils, ON, Canada, 1970; p. 274.

33. Launeau, P.; Archanjo, C.J.; Picard, D.; Arbaret, L.; Robin, P.Y. Two- and three-dimensional shape fabric analysis by the intercept method in grey levels. *Tectonophysics* **2010**, *492*, 230–239. [CrossRef]

34. Robin, P.-Y.F.; Jowett, E.C. Computerized density contouring and statistical evaluation of orientation data using counting circles and continuous weighting functions. *Tectonophysics* **1985**, *121*, 207–223. [CrossRef]

35. Krooks, A.; Kaasalainen, S.; Hakala, T.; Nevalainen, O. Correction of Intensity Incidence Angle Effect in Terrestrial Laser Scanning. *ISPRS Ann. Photogramm. Remote Sens. Spat. Inf. Sci.* **2013**, *II-5/W2*, 145–150. [CrossRef]

36. *Références Altimétriques Maritimes*; SHOM: Brest, France, 2017; ISBN 978-2-11-139469-8.

37. Dupuy, Q. *Cartographie Morpho-Sédimentaire Multisites et Evolution Spatio-Temporelles des Champs de Haploops sp. Mémoire de Recherche M2*; Université de Nantes: Nantes, France, 2018; 71p.

38. Cordier, C. *V/O Haliotis: Traitement des Données du Sonar Geoswath*; Rapport Interne Ifremer, Ifremer/DYNECO/EB/12-04/CC; IFREMER: Plouzané, France, 2012; 90p.

39. Augustin, J.; Edy, C.; Savoye, B.; Le Drezen, E. Sonar mosaic computation from multibeam echo sounder. In Proceedings of the OCEANS’94, Oceans Engineering for Today’s Technology and Tomorrow’s Preservation, Brest, France, 13–16 September 1994; Volume 432, pp. II/433–II/438.
40. Talley, L.D.; Pickard, G.L.; Emery, W.L.; Swift, J.H. Gravity Waves, Tides, and Coastal Oceanography. In *Descriptive Physical Oceanography*, 6th ed.; Swift, J.H., Ed.; Academic Press: Cambridge, MA, USA, 2011; Chapter 8.

<figure><img src="image_22.png" /></figure>

© 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).