###### Abstract
The CGD technique exhibits stable decomposition performance when applied to topographic LiDAR data in urban or forest areas where multiple return pulses are separate and have prominent peaks [23; 24]. However, it may not be suitable for asymmetric ABL waveforms with a peak left shift due to attenuation of the return pulse energy [25; 26], which can lead to point omissions for seafloor objects [18]. Different mathematical models for ABL waveform decomposition have been designed for different sensor and surveying environments. In most models, the ABL waveform is assumed to be decomposed into three distinguishable contributions: returns from the water surface, water column, and bottom (or seafloor). Returns from the water surface or bottom are generally approximated based on the Gaussian or Weibull model [14; 15; 26]. On the contrary, triangular [14], quadrilateral [15], and improved quadrilateral models [26], among others, have been used for modeling the returns from the water column considering the asymmetry due to the exponentially attenuated amplitude. Other studies have proposed mathematical models considering the water depth because it affects the shape of the returned signals. A depth-adaptive decomposition method was proposed by adopting two individual mathematical models according to the water depth [27]. One was an empirical function based on the calibration waveform for shallow water, and the other was an improved quadrilateral function with a second polynomial for deep water. A decomposition method for very shallow water (< 2 m) was also developed based on boxcar and exponential functions [28]. Waveform decomposition based on the separated functions for the three contributions can be efficient because the meaningful signals (i.e., the returns of the water surface and bottom for bathymetry) are simultaneously defined; therefore, bathymetric maps can be immediately generated. However, they are only applicable to typical waveforms in the ideal case (e.g., low-turbidity water without unexpected obstacles such as marine organisms or foam). In practice, the ABL waveform may not include returns from the bottom due to significant attenuation by the water column or include prominent returns from unexpected objects in the water, such as fish or seaweed. Multiple returns from the bottom due to seafloor unevenness may also be included. Furthermore, waveform decomposition techniques that rely on these fixed models encounter challenges when applied to further applications because of their limited ability to extract waveform features from each individual decomposed component. Hence, considering the varying characteristics of coastal and inland water, a flexible Gaussian decomposition method that does not fix the specified model for each contribution or limit the number of mathematical models may be desirable. A few methods employing a flexible number of Gaussian models for approximating returned ABL waveforms have been developed. One of the early studies on waveform data processing involved decomposing the ABL waveform into multiple Gaussian models and using the parameters of one of the extracted Gaussian distributions to estimate the seafloor roughness [9]. Another noteworthy study [25] proposed the Gaussian half-wavelength progressive decomposition (GHPD) method. The GHPD progressively decomposes a waveform using half-wavelength Gaussian functions based on the time sequence of received echo signals, aiming to reduce the problems of unreasonable decomposition and position shift of reflected pulse peaks caused by echo superposition. This approach is effective in decomposing ABL waveforms of various shapes and improves bathymetry performance, especially in shallow water. Other flexible Gaussian decomposition approaches were mainly investigated for application to topographic LiDAR data. Progressive waveform decomposition (PWD) was proposed to detect and fit distinct return components successively via peak detection [29].