Deep learning assisted exponential waveform decomposition for bathymetric LiDAR

The processing of bathymetric LiDAR waveforms is an important task, as it provides range and radiometric information to determine the precise location of water surface and bottom, and other characteristics like amplitude. The exponential waveform decomposition proved to be an effective algorithm for bathymetric LiDAR waveforms processing, however, it heavily relies on the high-quality initial estimates of the model parameters. This paper proposes to make use of deep learning to obtain the initial values directly from the input received waveforms without any hand-crafted features and prior-knowledges. Additionally, to provide training samples, we presents a method to create the synthetic bathymetric LiDAR waveforms by simulating of the backscatter cross function returned from water bodies. Two networks with different sensitivities of weak signals were trained by these synthetic waveforms, and used to estimate the initial values of the model parameters, a least square optimization follows up to obtain the final waveform decomposition result. This deep learning assisted exponential waveform decomposition method is applied to the real waveforms acquired by RIEGL VQ-840-G. The results show that estimations with the help of deep learning is less influenced by the intermediate peaks backscattered from objects and particles in water, producing a cleaner point cloud with less isolated points below water surface than the original exponential waveform decomposition. Moreover, the proposed sensitive DL-XDC is even able to detect some very weak bottom returns with low SNR.


Introduction
Airborne laser bathymetry (ALB) is an efficient method for underwater topography measurements.To deliver accurate 3D point clouds to the end users, the laser signals received by the detector need to be processed to obtain spatial and radiometric data.However, it's challenging to accurately interpret vast amounts of received waveforms across various water conditions.In an ideal water body, the laser pulse is expected to return one echo from the water surface and one echo from the water bottom.Due to the scattering of light by particles in the water, it is likely that a small portion of the pulse energy is reflected to the receptor when interacting with turbid water, resulting in one or multiple intermediate peaks apart from the echoes returned from the water surface and the water bottom.Additionally, the laser pulse is attenuated by absorption and scattering, which leads to only limited signals from the water bottom collected by the receiver.This makes it difficult to accurately distinguish the weak water bottom signal from background noise, especially in deep, turbid water.Therefore, a reliable waveform processing algorithm needs to be studied to retrieve precise range and radiometric features of the water body.
The received waveform is the convolution of the system response and the differential backscatter cross section (dBCS) of the target, which describes the physical properties of the target hit by the laser beam.The major task of our waveform analysis is to reconstruct the dBCS.Most methods that have been published w.r.t.ALB full waveform analysis, regarded the received waveform as a superposition of the signals of the water surface, the water column and the water bottom, in which the signals are described by mathematical functions.The signals of the water surface and the water bottom were represented as two Gaussian functions in the research of (Yang et al., 2022, Abady et al., 2014), while the authors in (Abdallah et al., 2013) fitted the water surface signal by a Gaussian function and the water bottom as a Weibull function.To include the system response into consideration, (Xing et al., 2019) used a transformation of the calibration waveform to fit the surface and bottom returns.The calibration waveform is the system response to a target with the Dirac-shaped backscatter cross section, like bare ground, and approximated by a smoothing spline function.
As for the fitting of the back-scattered water column signal, different shapes of the functions are used, such as the exponential function with a second-order polynomial used in (Xing et al., 2019), the quadrilateral function (Abady et al., 2014) and the triangle function (Abdallah et al., 2013).However, none of these methods considered or proposed a physical model to describe the laser beam backscattering within the water column.Waveforms collected by topography airborne LiDAR systems are normally decomposed into multiple Gaussian pulses under the assumption that the system response is a Gaussian, pulse and the convolution of two Gaussians yields a Gaussian as well.But this conventional Gaussian decomposition does not apply to the bathymetry waveforms with asymmetric shapes.Without consideration of the system response, the fitting methods that are merely applied on the received waveforms cannot fully capture the physical characteristics of the interaction of the laser beam and the water body.Additionally, a bias of the estimated water depth has been mentioned in (Abady et al., 2014) (Abdallah et al., 2013), which comes mainly from the approximated water column signals overlapping water surface and bottom returns, translating the peaks of the surface and bottom.
A deconvolution could be applied to remove the system response, the remaining waveform should then be the representation of the dBCS of the target.However, the deconvolution induces a large amount of noise as it amplifies the amplitudes of the high frequency noise-components, which causes the amplitude of these components that may contain meaningful information such as water surface and bottom to fall below the inherent noise of the system.The deconvolution in its bare form is unusable in practice.The authors in a comparative study (Wang et al., 2015) on the ALB waveforms processing state that the robust Richardson-Lucy deconvolution has a superior performance compared to the mathematical approximation of the received waveforms, as it follows the physical process of the signal.Their findings demonstrate the importance of including the system response in the ALB waveform processing.
To avoid the current shortcoming, an exponential decomposition with implicit deconvolution method has been proposed by Schwarz (Schwarz et al., 2017).The authors represented the dBCS of a water body by multiple exponential function segments.The received waveforms can then be modeled as a convolution result of the LiDAR instrument's system waveform and the dBCS of the illuminated target.The purpose of the waveform decomposition is to determine the values of the parameters of the dBCS model that can minimize the difference between the modelled waveform and the received waveform.As the solution to this non-convex problem is obtained by means of a modified Levenberg-Marquardt algorithm by imposing parameter constraints, the result is very sensitive to the choice of the initial values of the parameters.
In recent years, some research have investigated the potential of deep learning (DL) on LiDAR waveform processing.In the studies of (Lang et al., 2022) and (Fayad et al., 2021), the authors utilize a convolutional neural network (CNN) to retrieve canopy height from large footprint waveforms acquired by space-borne LiDAR.Their results have shown that the use of a CNN-based framework enables direct processing of waveforms, which avoids the necessity of computational manual algorithm calibration to find the best settings.Letard (Letard et al., 2023) designs a U-Net with additional attention layers to reconstruct the water surface, column and bottom components of the bathymetric lidar signal.Asmann (Asmann et al., 2021) regards the detection of principal peak components in waveforms as a classification problem, a 1D CNN was applied on waveforms and predicted the peak confidence of each elements in the waveform.This work is an extension of the exponential decomposition method proposed by Schwarz (Schwarz et al., 2017) and aims to utilize DL to estimate better than hand-crafted initial values of the dBCS model parameters for the exponential decomposition algorithm.The well-known ResNet-18 architecture is used and adapted to work on 1D waveforms.And the network is trained by synthetic waveforms.By interpreting the physical meaning of the parameters of the proposed dBCS model, we assign a reasonable range for each parameter and simulate the waveforms by randomly selecting values from the given range.In this way, a large amount of synthetic waveforms with ground truth parameter values can be created and fed into the designed network for the training.The estimations from the DL networks are provided as the initial values and then go through a nonlinear least square optimization to obtain the final values.To evaluate the performance, the proposed DL-assisted exponential decomposition is applied on real waveforms collected by a RIEGL VQ-840-G instrument to validate the performance.

Waveform modelling
The received waveform is considered to be a convolution of the system waveform and the dBCS of the target water body.The model of dBCS, denoted by   , consists of three contributions: water surface, water column, and water bottom.Based on the previous work (Schwarz et al., 2017), we use a Dirac-shaped pulse to describe the dBCS of water surface as it is considered a discrete scatter event.A truncated exponential function is used as the dBCS for water column, as both the laser pulse and the backscattered pulse attenuate exponentially when travelling in turbid water.Another truncated exponential function is used for the water bottom as it has proven to give a better fitting result than a Dirac-shaped dBCS.Thus, the dBCS of the water body, denoted by   can be written as the sum of the dBCS of water surface, water column, and bottom in Eq. ( 1).
Where () denotes Dirac's function, () 0 + denotes the Heavyside step function to indicate the effective time range of the dBCS of water column and water bottom.The Heavyside step function is given in Eq. ( 2).
The dBCS of water body (), is a function depending on 7 parameters (  ,   ,   ,  0 ,   ,   ,   ).The physical explanation of the parameters, which are also depicted in Fig. 1, are the following: The attenuation of the laser pulse travelling in water can be interpreted by the exponential function of the water column (the second item in Eq. ( 1)).The starting dBCS of the water column is   when the laser pulse hits the water surface, and drops to the value of   / after travelling   in the water.In turbid water, the energy decreases quickly due to strong scattering by dense submerged particles in the water, hence the value of   is very likely to be low.While in clear water, the laser pulse tends to travel for a long distance without so much energy attenuation, hence the value of   is very likely to be high.
Although it could be expected that the response of the water bottom would be a Dirac-shaped function, the dBCS of the water bottom of this work uses an exponential function for a better approximation of the received waveform.Therefore, the physical meaning of the water bottom PRD (  ) and decay time (  ) cannot be interpreted as we did for the water column, these two parameters are merely used for mathematical approximation.
The returned energy of the three target components (surface, column, and bottom) can then be calculated by the integral of their dBCS functions over the according time range.As plotted in Fig. 1, the two shaded areas under the two curves correspond to the returned energy of water column and water bottom, respectively.The mathematical equations are given in Eq. (3), Eq. ( 4) and Eq. ( 5), where   ,   ,   denote returned energy from the water surface, the water column, the water bottom, respectively.  ,   ,   are in units of 1. Figure 1.An example model of the dBCS of a target water body that contains waveforms returns of water surface, column and bottom.
The received signal   can be written as a convolution with   the dBCS and the system waveform ℎ  , as shown in Eq. ( 6).
The aim of the exponential waveform decomposition (XDC) is to iteratively compute the 7 parameters of the dBCS function in Eq. ( 1), so that the modelled waveform, which is a convolution (computed by Eq. ( 6)) between the given system waveform and the estimated dBCS, is best-fitting the received waveform.In the processing of the original XDC (Schwarz et al., 2017), the initial values of the dBCS model parameters need to be firstly estimated by a few hand-crafted features, based on plausible reasoning, then a least square method is applied to progressively determine the final values.The result of the least square optimization heavily depends on the proximity of the initial values to the solution to avoid being trapped into a suboptimal result.In this paper we propose to replace the hand-crafted method by DL to facilitate the initial estimates of the dBCS parameters.We use the DL network to perform the first step in the XDC, which outputs the parameter values directly from the input waveforms, then go through the same least square procedure to obtain the optimized values of the dBCS.

Waveform simulation
To train a proper deep learning network, a large amount of accurately labelled training data is required.We used simulated waveforms as training data.As the dBCS model of water body is driven by the physical characteristics of the target water area, the dBCS of water body can be simulated by setting a reasonable range for the parameters and convolved with the system waveform of the specific scanner device to obtain synthetic echo waveforms.The value ranges are presented in Tab. 1 and detailed as below.

•
The total length of the waveform is set to be 512 ns.The start time ( 0 ) is related to the instrument's height above the water surface and the beginning of the waveform recording.Thus a wide range is chosen for  0 as the water surface echo could appear at any location of the waveform.

•
The time difference (  ) between the surface and the bottom return is related to water depth.The possible value of   ranges from 1.5 ns to 250 ns, which corresponds to the water depth varying from 0.165 m to 27.5 m.The half speed of light travelling in water (0.11 m/ns) is used to convert time to range.

•
The decay time (  ) of water column dBCS relates to the water turbidity.It tends to be high for a clear water body, whereas it is low when water is turbid.Based on Eq.
(1), the starting dBCS of the water column is   , and drops to   / 1/  after travelling 1 ns in the water, thus the drop rate is 1/ 1/  .We set the range of   as 0.8 ns to 100 ns.That is to say, in an extremely turbid water, the dBCS drops to 28% while the dBCS can keep 99% in an extremely clear water after the laser pulse travelling 1 ns in the water.

•
The decay time (  ) of water bottom dBCS should be a relatively small value, as the water bottom is solid target and the laser pulse should be mostly back-scattered.It uses the same minimum values as the   used, and maximum values is set empirically.

•
Instead of interpreting the range values for PRD, we computed the PRD of the water column and the water bottom by its returned energy and decay time based on Eq. ( 4) and Eq. ( 5), because the returned energy of the target is easier to interpret than the PRD.

•
The maximum value is set as 0.88 in units of 1, for the returned energy of the water surface (  ), the water column (  ) and the water bottom (  ).Two different minimum values are used for the returned energy of the water bottom to train two networks with the different sensitivities of the weak signals.
The returned energy of water column and water bottom are negatively correlated.When the energy returned from the water column is high, for instance, in a turbid water body, the returned energy from the water bottom should be low.In a clear water body, the water bottom energy tends to be high but results in a low energy returned from water column.Therefore, we added a constrain in the simulation as shown in Eq. ( 7) to generate more realistic waveforms.We simulated 6 types of waveforms by varying the returned energy of the water surface, column and bottom.
For each type of waveform, parameter values are randomly selected from the given range with the constrain in Eq. ( 7) to create various dBCS functions, the synthetic waveforms are then generated by convolution of the given system waveform and the simulated dBCS functions.
The real waveform received from the LiDAR sensor is very unlikely to be as smooth as the simulated waveform, thus we considered three types of noise to complement the simulation, which are random Gaussian noise, one or two spike-like outliers occurring at arbitrary location, and a false-alarm signal with a very small energy before the water bottom return.This false-alarm signal was generated by convolving the system waveform and a small Dirac-shape pulse to simulate the peaks that are caused by backscattering at particles.To improve the robustness of the trained network, the three types of noise were added during the data augmentation of the network training.

Deep learning network
The ResNet-18 architecture (He et al., 2016) is adapted in this work.All the two-dimensional operations are replaced with one-dimensional convolution, pooling and batch normalization layers, in order to deal with 1-D input waveforms.The overall architecture of the network is presented in Fig. 2(a).The ResNet-18 comprises 8 residual blocks, as depicted in Fig. 2(b), each block consists of two consecutive 3x3 convolutions, each followed by a batch normalization and a rectified linear unit (ReLU) to introduce non-linearity.In our work, the first convolution in the 1st, 3rd, 5th and 7th block is performed with a stride size of 2. This down-sampling operation gradually increases the receptive field of the input feature vector, so that the parameters that need a large neighbourhood to interpret, like  0 , could be well learned through the features represented by deeper layers.
After the last residual block, a global average pooling layer is added, to average all sample index in each feature dimension to a scalar value, and result in a 1-dimensional feature vector.At last, a fully connected layer is applied to output the 7 target parameters' values.Note that the 7 target parameters are the start time ( 0 ), the time difference (  ), the decay time of water column (  ), the decay time of water bottom (  ) , the energy of water surface (  ), the energy of water column (  ) and the energy of water bottom (  ).
The water column PRD (  ) and the water bottom PRD (   ) are calculated based on Eq. ( 4) and Eq. ( 5).We simulated 70,000 waveforms of each type and use the resulting 420,000 waveforms to train the networks.Two networks with different minimum values for the water bottom energy are trained.The waveforms that are simulated with the relatively higher energy threshold of the water bottom return are used to train the network that would be more robust against noise signals, named robust-DL network in the following.Whereas, the network trained by the waveforms with a relatively lower energy threshold of water bottom is expected to be more sensitive to the weak bottom returns.

Results
The performance of the proposed method was tested on real waveforms acquired by RIEGL VQ-840-G in the bay area of Toulon, a Mediterranean coastal city in France.Due to the lack of ground truth values, we compare the results against the results of the original implementation of exponential decomposition, which is denoted as Ori-XDC.In the Ori-XDC, the initial values of surface and bottom returns were computed by evaluating the isolation, prominence and amplitude of each sample point.The details can be seen in (Schwarz et al., 2019).DL-XDC denotes the exponential decomposition where initial values were estimated by a deep learning network.Based on which pretrained network has been used for initial values estimation, we refer to our methods as sensitive DL-XDC and robust DL-XDC.Note that the same modified un-constrained least square optimization method is applied after estimating the initial values in the above-mentioned methods.The results are evaluated by visually inspecting the generated point clouds in Section 3.1, the waveform decomposition results in Section 3.2 comparing the modelled waveforms and the received waveforms and the computation time in Section 3.3.

Generated point clouds
Fig. 3 presents the point clouds generated by each decomposition method.The length of the profiles in Fig. 3 is about 162 m and the deepest area is about 29 m deep.The Ori-XDC tends to generate more points within water column.Compared with robust DL-XDC, more bottom points are detected by the sensitive DL-XDC method, however, some points under the water surface are also returned due to its high sensitivity to weak signal returns.For better evaluation of the bottom detection, we manually identified the bottom points and had the selected points counted in RiPROCESS (RIEGL, 2018).The resulting bottom point count by Ori-XDC, robust DL-XDC sensitive DL-XDC is 458, 560 and 630, respectively.This suggests that an improvement of bottom detection is achieved by both DL-XDC methods, meanwhile resulting also in less points within the water column.

Waveform decomposition results
To intuitively examine the effect of DL, we selected two received waveforms that are difficult to interpret by handcrafted features such as peak isolation and prominence, and depicted their modelled waveforms to illustrate the decomposition performance.The modelled waveforms are the result of the convolution of the system waveform and the target dBCS function, in which the coefficients are solved by Ori-XDC and the proposed DL-XDC methods.The intermediate modelled waveform before the least square optimization procedure are presented as well, in order to evaluate the quality of the initial values provided by Ori-XDC and DL methods.
Since the laser pulse may be scattered by objects within the water column, an intermediate prominent peak, sometimes with a even stronger returned energy than the water bottom echo, may be observed in the received waveform, as for example in Fig. 4. The water bottom echo is initially detected at the index of 472, 621 and 620 by Ori-XDC in Fig. 4(b), robust-DL in Fig. 4(c) and sensitive-DL in Fig. 4(d), respectively, After the iteratively optimization, the water bottom echo is determined at the index of 472, 622 and 622 by Ori-XDC in Fig. 4(e), robust DL-XDC in Fig. 4(f) and sensitive .By hand-crafted features, Ori-XDC identified the intermediate peak as the bottom return because it is the second most prominent peak detected in the waveform, which results in many scattering points (Fig. 3(a)) below the water surface.Both DL methods tend to regard the last weak signal as the bottom return, thus yielding less scatter points in the water column.
In Fig. 5, the water bottom return is very weak and and nearly obscured by noise.For the initial estimates, the bottom echo is detected at the index of 429, 683 by Ori-XDC in Fig. 5(b) and sensitive-DL in Fig. 5(d), respectively, while no bottom echo is detected by robust-DL in Fig. 5(c).After the optimization, the bottom echo location is finetuned and detected at the index of 421 and 688 by Ori-XDC in Fig. 5(e) and sensitive DL-XDC in Fig. 5(g), respectively.As the optimization is only performed on the input echo parameters, no new echo would be created during the optimization procedure.Thus, only the surface echo parameters are fine-tuned for the robust-DL.The sensitive-DL network was trained to be more responsive to weak signals and successfully detects the very weak water bottom signal from such a low signal-noise-ratio (SNR) waveform.As for Ori-XDC, an intermediate peak is considered as the water bottom echo like in the previous example shown in Fig. 4, The waveforms modelled by the initial estimates of the DL models are plotted in the (c) and (d) in Fig. 4 and Fig. 5.The modelled waveforms fit the received waveforms already quite well before the optimization procedure proving that the trained DL networks provide more accurate estimates of the initial values.
Considering the test area is a relatively clear water body, the water column is expected to be clean rather than occupied by points.Additionally the water bottom points generated by the DL methods are in agreement with the neighbouring bottom points.Thus we conclude that the DL methods improve the waveform decomposition performance and result in more accurate water bottom points detection.Both DL-networks show the ability to better distinguish intermediate signals back-scattered by objects or particles and the bottom returns.In addition, the sensitive DL network is able to detect even weak bottom returns exhibiting very low SNRs.

Computation time
Least square optimization iteratively minimize the difference between the received waveform and the modelled waveform to determine the best parameter values of the dBCS function.Better initial values result in fewer optimization steps and thus likely less processing time.Tab. 2 presents the processing time for each method for the record containing 670309 waveforms.Runtime denotes the overall processing time from reading the input data to writing the point cloud, whereas full waveform analysis time is only the time spent on waveform decomposition itself.
For Ori-XDC, the waveform decomposition time is 22m47s, which includes the initial values estimation and the optimization procedure.
For DL-XDC, the waveform decomposition comprises two steps, the initial values estimation by the trained DL networks and the follow-up optimization step implemented by XDC algorithm.Note that the deep learning inference runs on a GPU (Quadro M2200).

Conclusion
In this work we designed and trained a deep learning network to provide initial values for the parameters of the dBCS used for full waveform analysis in the context of bathymetric laser scanning.The proposed DL-assisted waveform decomposition method is evaluated on the ALB dataset collected in the ocean area.The improvement of these initial values compared to hand-crafted initial values can be confirmed by the generated point clouds, modelled waveforms by the estimated parameters, and the computation time.Compared to the Ori-XDC, the resulting point cloud with the assistant of DL has less isolated points in the water volume but more water bottom points.By visually inspecting the waveform decomposition results, the DL-XDC methods are less influenced by the intermediate peaks in the received waveform and more likely to detect the real water bottom return.Two DL-XDC networks are offered based on different settings for the SNR: more bottom signals can be detected by the sensitive DL-XDC in deep water areas, but generating relatively more scattering points within the water column than the robust DL-XDC.In the future work, the performance of the proposed DL-XDC will be evaluated on the ALB waveforms collected in different water bodies, such as lakes and rivers.bottom-detection algorithm for LiDAR bathymetry of very shallow waters.ISPRS Journal of Photogrammetry and Remote Sensing, 150, 1-10.Schwarz, R., Pfeifer, N., Pfennigbauer, M., Ullrich, A., 2017.Exponential Decomposition with Implicit Deconvolution of Lidar Backscatter from the Water Column.PFG -Journal of Photogrammetry, Remote Sensing and Geoinformation Science, 85, 159-167. Wang, C., Li, Q., Liu, Y., Wu, G., Liu, P., Ding, X., 2015.A comparison of waveform processing algorithms for single-wavelength LiDAR bathymetry.ISPRS Journal of Photogrammetry and Remote Sensing, 101, 22-35. Xing, Wang, Xu, Lin, Li, Jiao, Zhang, Liu, 2019 DL network (a) 1D ResNet-18 used for waveform parameters estimation, (b) Residual Block sensitive DL-XDC Fig. 3. Near shore sea profiles of the point clouds.The blue point clouds are the echoes from the water surface and the red point clouds are the echoes from the water bottom.(a) was generated by Ori-XDC, (b) was generated by robust DL-XDC and (c) was generated by sensitive DL-XDC.
example of the received waveform with the echo returned from submerged particles or objects in water.(a) received waveform, the subset in the green rectangle are further depicted with additional modelled waveform in (b)-(g) to have a close look at the fitting results.(b)-(d) are the modelled waveforms (before least square optimization) constructed by the initial parameters values estimated by Ori-XDC, the predicted parameter values directly from robust DL and sensitive DL, respectively.(e)-(g) are the modelled waveforms after least square optimization using the initial values provided by the different methods.

Fig. 5 .
Fig. 5.An example of the received waveforms with the weak signal of the water bottom.(a) received waveform, the subset in the green are further depicted with additional modelled waveform in (b)-(g) to have a close look at the fitting results.(b)-(d) are the modelled waveforms (before least square optimization) constructed by the initial parameters values estimated by Ori-XDC, the predicted parameter values directly from robust DL and sensitive DL, respectively.(e)-(g) are the modelled waveforms after least square optimization using the initial values provided by the different methods.
Thus, the total waveform decomposition for DL methods is the sum of initial values estimation time by DL and optimization time by XDC.The waveform decomposition of robust DL-XDC and sensitive DL-XDC is 18m57s and 18m38s, respectively.Overall, the waveform decomposition time required by the DL-XDC is slightly less than the Ori-XDC.
. A Depth-Adaptive Waveform Decomposition Method for Airborne LiDAR Bathymetry.Sensors 19, 5065.Yang, F., Qi, C., Su, D., Ding, S., He, Y., Ma, Y., 2022.An airborne LiDAR bathymetric waveform decomposition method in very shallow water: A case study around Yuanzhi Island in the South China Sea.International Journal of Applied Earth Observation and Geoinformation, 109, 102788.