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

: 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 ﬁrst derivative of a wide Gaussian ﬁlter serves to delineate the full-waveform range limits. Because of turbidity current heterogeneity and complexity of multiple layers radiative transfer model, a drastic simpliﬁcation 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.


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][2][3][4][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].
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.  [23]); red +1 and +6 are sample locations of data used for presentation of the method.
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 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.

Full-Waveform Data Registration and Ray Tracing
Like many authors who have studied the diffusion of light through the thickness of vegetation [24][25][26][27][28][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].

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 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): Remote Sens. 2019, 11, 117 6 of 23 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 Figure 2b. 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 Figure 1f): (a) Initial FWF relative intensity DN versus range R presented between Rb and Re 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.
Following the Snell-Descartes law, as in [13] for example, light speed slows down from air ca to water cw and the incident angle of a laser beam in air aa rotates to aw, according to their respective refractive indices na and nw. Since index na is equal to 1 the 532 nm laser beam in water slows down proportionally to 1/nw and the incident angle in water becomes arcsin(sin aa/nw). 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  Figure 1f): (a) Initial FWF relative intensity DN versus range R presented between Rb and Re 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.
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 (α 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 α w and the cosine of incident angle in air α a . Therefore, at a 1-m-wide pixel location with an incident angle α a in air with a 28 • FOV, the full vertical range water index correction becomes: Much higher incident angles would require wave slope corrections as shown by [35], which also justifies our choice of 28 • FOV.
Remote Sens. 2018, 10, x FOR PEER REVIEW 7 of 23 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 (αi) and its range in water (Rw), xa and xw respectively in Figure 3, is equal to the ratio between the cosine of incident angle in water αw and the cosine of incident angle in air αa. Therefore, at a 1-m-wide pixel location with an incident angle αa in air with a 28° FOV, the full vertical range water index correction becomes: Much higher incident angles would require wave slope corrections as shown by [35], which also justifies our choice of 28° FOV.  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 m −1 as shown in Figure 3b. Moreover the 3 × 3 pixel median used to reduce the noise retaining only 1 best out of 9 neighboring pixels the effective resolution is 9 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. 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. 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. 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 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 T 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 T and the last echo must be determined for each study. It is 2.  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 T 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 T can be used to infer the range of the last unshaped echo with the addition of an offset determined by ground control points.

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.
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.  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 T 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 T can be used to infer the range of the last unshaped echo with the addition of an offset determined by ground control points.

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.
Remote Sens. 2018, 10, x FOR PEER REVIEW 9 of 23 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 T 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 T 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, Figures 2a,d and 5b. Its full determination is part of Section 3 using GPS ground control points. 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 T 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 T can be used to infer the range of the last unshaped echo with the addition of an offset determined by ground control points.

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.
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.  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. 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.

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 T 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 T bottom in Figure 8c.
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   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. 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

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.

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.

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 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 T , (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 T 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 T bottom forming a less smooth profile. 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 11 a and 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 T , (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 T 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 T bottom forming a less smooth profile.  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 T 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.

GPS-derived z elevations, (1) in
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 T 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 T ) −1 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 T ) −1 data in the following work.
Remote Sens. 2018, 10, x FOR PEER REVIEW 14 of 23 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 T 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 T 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 T ) −1 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 T ) −1 data in the following work.  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. 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 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 lowpasses used to minimize them. 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 T bottom with the second method of bathymetry mapping.

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 6m 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 T 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 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. 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 T bottom with the second method of bathymetry mapping.

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 6m 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 T 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 T , (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.
Remote Sens. 2018, 10, x FOR PEER REVIEW 16 of 23 water index corrected NCFWF T , (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.  water index corrected NCFWF T , (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.  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.

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 2 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. 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 (Figure 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.  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 (Figure 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.  Despite a stronger noise level than echoes obtained with the dddNCFWF first method, the NCFWF T 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. Despite a stronger noise level than echoes obtained with the dddNCFWF first method, the NCFWF T 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.

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

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 16c.
The 2.8 m offset between last echo and the bottom of the FWF damping function given by NCFWF T 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 T 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 T 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. 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). 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 2 . 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.