## 3 Methods  
The green laser beam used for ABL traverses the water surface and propagates in the water column until it reaches the bottom (Figure 2). The amplitude is rapidly attenuated by absorption, scattering, and refraction. If the beam is not entirely attenuated while traveling through the air-water interface, the backscatter signal reaches back to the receiver,  
Figure 1: Test sites of Seahawk: (**a**) locations of the test sites; (**b**) test site 1 on eastern coast of Donghae, Korea; (**c**) test site 2 on southeastern coast of Jeju Island, Korea.  
revealing the water surface and bottom positions through peaks in the waveform [35]. A typical ABL waveform is composed of the water surface return, water column backscatter return, bottom return, and noise [10]. Each component is separated through waveform decomposition, and points from the decomposed components are registered.  
A LiDAR waveform is a convolution between a transmitted laser pulse and a surface scattering function, both of which are generally considered to follow a Gaussian model. Waveforms received through multipaths can be represented as Gaussian mixture models [24], as in Equation (1), where \(n\), \(A_{i}\), \(\mu_{i}\), and \(\sigma_{i}\) denote the number of Gaussian models, amplitude, center position, and standard deviation of the ith Gaussian model, respectively. Because the parameters of Gaussian mixture models have physical meanings related to the position, intensity, and properties of the signal, the characteristics of individual signals can be analyzed through the parameters. To approximate the waveform to the nonlinear model, iteration optimization, such as the Levenberg-Marquardt technique, is performed. In this process, initial values for each parameter (\(n\), \(A_{i}\), \(\mu_{i}\), and \(\sigma_{i}\)) are required. In particular, it is important to select initial values for the iterative operation to finally converge to a global rather than a local solution.  
\[\hat{y}=\sum_{i=1}^{n}A_{i}exp\left(-\frac{\left(t-\mu_{i}\right)^{2}}{2\sigma _{i}^{2}}\right) \tag{1}\]  
Generally, the Gaussian decomposition technique approximates a waveform to the Gaussian mixture models using peak parameters (amplitude, center, and width) extracted via peak detection as initial values. This CGD relies on peak detection results and is suitable for waveforms in which individual return pulses are detected as separate peaks. However, the ABL waveform is continuously scattered during underwater signal propagation, and multiple Gaussian return signals are densely overlapped, resulting in a signal tilted to the left. It is challenging to separate the overlapped underwater backscattering components by CGD, which consequently imposes limitations on precisely decomposing the water surface and bottom returns.  
In this study, a progressive Gaussian decomposition (PGD) method is proposed to effectively decompose the ABL waveform data. The proposed approach gradually decomposes a waveform by estimating the overlapping peaks that are undetected during the initial peak detection (Figure 3). The PGD method uses the detected peaks as the first initial  
Figure 2: Data generation in airborne bathymetric LiDAR system.  
value but is not bound by them. It repeats multiple Gaussian curve fitting, progressively estimating potential peaks (PPs) until the desired criteria are achieved. Denoising and signal range selection are performed as preprocessing in the proposed approach. The detailed process is as follows.