Gaussian Half-Wavelength Progressive Decomposition Method for Waveform Processing of Airborne Laser Bathymetry

Abstract

In an airborne laser bathymetry system, the full-waveform echo signal is usually recorded by discrete sampling. The accuracy of signal recognition and the amount of effective information that can be extracted by conventional methods are limited. To improve the validity and reliability of airborne laser bathymetry data and to extract more information to better understand the water reflection characteristics, we select the effective portion of the original waveform for further research, suppress random noise, and decompose the selected portion progressively using the half-wavelength Gaussian function with the time sequence of the received echo signals. After parameter optimization, a reasonable and effective reflection component selection mechanism is established to obtain accurate parameters for the reflected components. The processing strategy proposed in this paper reduces the problems of unreasonable decomposition and the reflected pulse peak-position shift caused by echo waveform superposition and offers good precision for waveform decomposition and peak detection. In another experiment, the regional processing result shows an obvious improvement in the shallow water area, and the bottom point cloud is as accurate as the intelligent waveform digitizer (IWD), a subsystem of airborne laser terrain mapping (ALTM). These findings confirm that the proposed method has high potential for application.

1. Introduction

An airborne LiDAR bathymetry (ALB) system is a type of full-waveform scanning bathymetry system based on an aerial platform. This technology has been increasingly used in the field of bathymetry because of its features of stability, accuracy, and efficiency [1]. The full-waveform data processing technique is the basis of ALB systems, such as the Aquarius and CZMIL systems made by Teledyne Optech [2,3], the VQ–820G and VQ–880G series produced by Riegl [4,5], the LADS Mk series produced by Fugro [6], and the Hawk Eye series and Chiroptera II produced by Leica [7,8]. In the process of laser propagation, energy and direction are usually affected by the condition of the medium because laser propagation is often accompanied by reflection, refraction, scattering, and absorption, especially at the interface between air and water. The reflected signal received by an ALB system, which usually samples discretely, exhibits a complex waveform shape because variations in the water-body bottom produce different backscattering characteristics [9].

In full-waveform data processing, improving the signal-to-noise ratio and accurately identifying surface and underwater temporal positions are prerequisite problems of ALB systems. Enhancing the ability to identify the reflected signals of water surfaces and bottoms has great significance for the entire system, especially in shallow or turbid water [6,10]. Billard proposed an algorithm for detecting the submarine reflection signal in an echo signal, in this algorithm, a high-pass filter signal is used to remove the low-frequency part of the echo signals, and two types of methods for identifying the high-frequency pulse have been proposed [11]. Henry Wong and Andreas Antoniou proposed a digital smoothing filtering method to eliminate the signal noise while effectively retaining the echo signal [12], and they also proposed a method of using a low-pass digital differential filter to extract sea bottom reflection signals [13]. In their research, they considered that if the filter could select the cut-off frequency combined with the physical and chemical properties of the water bottom, which can influence the shape of echo waveforms, then the surface, column and bottom of the water could be distinguished to some extent from the complex echo waveform. This method has also been experimented with by other researchers and has been verified and applied in their related experiments [14]. Taubin proposed a $\lambda /\mu$ filtering algorithm that defines a low-pass filter based on the Gaussian function and can effectively avoid the data filtering edge contraction to retain the edge information of the original data [15]. The effective part of the waveform signal after noise smoothing contains two reflections, and the energy changes when the laser penetrates the water. The main purpose of waveform processing is to detect the temporal positions of reflections, to identify the echo type and to extract the relevant information. Because the detection environment is always complicated due to potential variations in the properties of the benthic layer, the water transparency, and the depth, all the requirements cannot be met by using a perfect echo-data processing algorithm [16]. Conventional waveform processing methods can be categorized into three types:

• Echo detection uses the shape of the echo waveform to identify the temporal location of the energy mutation that is thought to be the occurrence time of the reflection. Conventional methods of echo detection include the maximum peak (MP), zero-crossing, and averaged square difference function (ASDF) methods [17].
• Deconvolution is usually applied for image or signal restoration. Jutzi and Stilla first confirmed that this method can effectively extract a target with a height greater than 15 cm [9]. The Gaussian decomposition method and the deconvolution method have been compared, and the deconvolution technique was found to obtain more peaks than the decomposition method [18]. Wang compared the results of several conventional methods in single-band laser data processing and found that the Richardson–Lucy deconvolution method has a higher detection rate and lower error. The disadvantage of deconvolution is that its anti-noise ability is always weak and can easily cause the misjudgment of the echo temporal location. Moreover, this method is also prone to ringing effects that adversely affect the results of data processing [19].
• Mathematical approximation. Generally, ALB system samples the reflected signal at a certain frequency. Hofton proposed that the echo waveform should be understood as the superposition of several Gaussian components [20]. Several Gaussian decomposition methods have been widely used, such as layered Gaussian function waveform fitting based on nonlinear least squares [21], the Gauss–Newton method [22], and the EM algorithm [23]. Both Zwally and Wagner suggested that the results obtained by using Gaussian decomposition are more consistent with the needs of multidisciplinary applications [24,25]. Moreover, the multiscale wavelet analysis method was used in waveform decomposition for light detection and ranging waveform characterization, and the results shows the consistency with the GLA14 product [26].

Most ALB system manufacturers do not currently disclose their processing methods for full-waveform data because of commercial confidentiality. The following problems are observed in waveform data processing for applications of ALB technology.

• ALB systems collect intensity data from the amplitude of received signals by discrete sampling at a certain frequency; however, the full-waveform data contain a considerable amount of redundancy. Operationally, the effective part of the echo must first be locked to reduce misjudgments and to improve computational efficiency. Additionally, the noise contained in the waveform data has adverse effects on the accurate determination of the reflection time; thus, smoothing filtering should be performed according to the echo characteristics before analyzing the original waveform data.
• Determining the temporal positions of abrupt changes in the echo energy by directly using the discrete sampling points is difficult. The simple interpolation results are always in error with respect to the actual reflection time [20]. A smooth curve fitted from discrete full-waveform data would facilitate determining the temporal position of abrupt events and the propagation time between the different reflected waveforms.
• An echo signal is not a regular Gaussian function, and its shape is affected by the attenuation of the medium, which often shows a trailing characteristic. When the distance between reflected objects is small or the water depth is shallow, the received waveforms may overlap. This situation shifts the apparent bottom reflection, and the corresponding zero crossing does not represent the true bottom-peak position. The worst case is that the bottom-reflected signal is completely embedded in the echo signal, which will increase the complexity of water depth estimations [18]. Therefore, a reasonable component selection mechanism must be adopted.

To address the problems of full-waveform data processing, a full-waveform simulation method for Gaussian half-wavelength progressive decomposition (GHPD) is presented in this paper. We used this method with full-waveform data received by the ALB system Aquarius as an example and overcame the problems of low waveform-simulation fitting accuracy, the lack of clarity in the process of producing the signal, and peak-position offsets caused by overlaps. We compared the effects of the GHPD method to those of the intelligent waveform digitizer (IWD), a digitization and recording system produced by Teledyne Optech, and analyzed and summarized the characteristics of the GHPD method in full-waveform data processing to provide a necessary reference for further research and improvements.

2. Materials and the Methods

2.1. Condition of the Experimental Area

The Aquarius system was produced in 2011 by Teledyne Optech, and it was based on the Gemini system, which was suitable only for shallow water areas. The effective design detection range is 15 m, and the Aquarius system uses linear scanning to cover the target area, with a scanning angle of ±20°. To verify the performance of the ALB system and its applicability in the China Sea, the First Institute of Oceanography of the State Oceanic Administration (SOA) conducted depth measurement experiments together with Teledyne Optech in the South China Sea at the end of 2012.

We selected partial airborne laser bathymetry data from two different types of bottom material as the experimental data to represent two classical seafloor types (Figure 1). Experimental area A is located in the northwest portion of the island. The area is shallow and gentle, and the bottom is fine sand. The average depth within the perimeter of the island is approximately 6.5 m. The reflected signals are also received in areas with depths greater than 9 m. Experimental area B is located in the sea reef site; coral reefs and coral reef debris are the main material on the bottom, so the seabed topography is complex and changeable. The farthest offshore distance of Experimental area B is 1240 m, and the maximum depth is approximately 20 m. The transparency of both experimental areas is approximately, the Secchi depth is 8 m, and the sea conditions were stable on the day of the experiment.

2.2. The Full-Waveform and the GHPD Process

In the case of a simple scatterer, if complex noise interference and signal attenuation are not considered, then the echo signal is a superposition of the signals generated from the different reflection cross-sections and Gaussian white noise. The reflected pulse energy that is returned approximately obeys the Gaussian distribution, and the full-waveform data can be expressed as follows:

$$x(t)=n(t)+{\displaystyle \sum _{i=1}^m{f}_i(t)}=n(t)+{\displaystyle \sum _{i=1}^m{A}_i{e}^{-{(t-{\mu }_i)}^2/2{\sigma }_i^2}}, n(t)~N(0,{\sigma }^2),$$

(1)

where $f_i(t)$ is the time response function of each component in the full-waveform data and $n(t)$ is the Gaussian white noise. The three parameters of the Gauss function $(A_i,{\mu }_i,{\sigma }_i)$ are unknown, and they represent the amplitude, temporal position and waveform width of the reflected waveform of each reflection section, respectively. Obviously, the essence of Gaussian fitting is to determine the values of these three parameters. In this paper, the proposed GHPD method works by observing the shape of the echo waveform. Figure 2 shows the main process of the GHPD method.

2.3. Data Pretreatment

Data pretreatment in this paper refers to the adjustment and preparation of the full-waveform data to be processed. To reduce data redundancy, determine the effective segmentation of the waveform and reduce the effects of noise in the process of laser propagation, two steps are performed: selection of the effective part and suppression of the original waveform.

2.3.1. Selection of the Effective Part

The effective part of the echo waveform is the section of amplitude change within the entire waveform caused by reflection. Concretely, intensity sequences of the laser signal exhibit acute and obvious fluctuations in the entire waveform.

(1) Echo waveform segmentation

The term $data(t)$ represents the original full-waveform signal sequence, and it can be divided into three parts based on the time sequences of the received pulse signal: $x_1(t)$, $x_2(t)$, and $x_3(t)$:

$$\{\begin{array}{l}{x}_1(t_1)=data(t_1), t_1<t_a\\ x_2(t_2)=data(t_2), t_a\le t_2\le t_b\\ x_3(t_3)=data(t_3), t_3>t_b\end{array},$$

(2)

where $t_a$ and $t_b$ are the starting and ending times, respectively, of the effective part $x_2(t_2)$.

(2) Selection principle of the waveform effective part

The shape of the echo waveform shows that the tail part of the signal is gradually reduced by the influence of the laser pulse energy, while the noise intensity is similar to that of the former part (Figure 3). Therefore, the statistical characteristics of the effective part of the waveform are obviously different from those of $x_1(t)$ and $x_3(t)$. We use three times the standard deviations of $x_1(t)$ and $x_3(t)$ as respective thresholds and constantly add sampling points from both edges of the waveform to the effective part to update the amplitude sequences in each side, The amplitude sequences in the same side include two sequences with different ends, that one stops at the temporal position corresponding to the local maximum, and the other stops at the temporal position of the corresponding to the local minimum. Before the new standard deviations of the sequences with the local maximum end are greater than the thresholds, the parts of the waveform sequences mainly contain the effects of reflected pulses and the interference of background noise. Therefore, we truncate the waveform data, and define the effective part between the $t_a$ and $t_b$. $t_a$ and $t_b$ are the temporal positions of the respective located minima.

(3) Method for selecting the effective part of the signal

To truncate the effective part at the local minimum, we use the first derivative to calculate the local extreme of the waveform and obtain the amplitude value to build the two amplitude sequences from the beginning ($t=t_0$) and the end ($t=t_{end}$) of the waveform to the effective part.

Local maximum sequence: $x_u(t)$, $t=t_m^u$, where $t_m^u$ is the corresponding temporal position sequence and $x_d(t)$ represents the total number of the maximum sequence.

Local minimum sequence: $x_d(t)$, $t=t_k^d$, where $t_k^d$ is the corresponding temporal position sequence and $t_a=t_{dj}$ represents the total number of the minimum sequence.

The interval between the maximum and the minimum is

$${\epsilon }_{m-k}=-{\epsilon }_{k-m}=t_m^u-t_k^d.$$

(3)

Therefore, the standard deviation of the extended sequence in different starting directions is given by

$$\{\begin{array}{l}{\sigma }_u(m)=\sqrt{\frac{{\displaystyle \sum _{t=P}^{t_m^u}{\left[x(t)-\overline{x}\left(t_m^u\right)\right]}^2}}{\left|t_m^u-t_0\right|}}\\ {\sigma }_d(k)=\sqrt{\frac{{\displaystyle \sum _{t=P}^{t_k^d}{\left[x(t)-\overline{x}(t_k^d)\right]}^2}}{\left|t_k^d-t_0\right|}}\end{array},$$

(4)

where

$$\overline{x}(t_m)=\frac{1}{\left|t_m-P\right|}{\displaystyle \sum _{t=P}^{t_m}x(t)}.$$

(5)

Left to right ($P=t_0$): $m=i$, $i=1,2,\dots N_u$ and $k=j$, $j=1,2,\dots N_d$; Right to left ($P=t_{\mathrm{end}}$): $m=N_u-i$, $i=0,1,2,\dots N_u-1$, and $k=N_d-j$, $j=1,2,\dots N_d-1$.

In the process of updating the sequence, the local amplitude values would be added into the sequence until it reaches the corresponding temporal position of the next extremal value. If the standard deviation of the new sequence with local maximum end is three times greater than the standard deviation of the sequence with minimum end, the effective part of the reflected signal is considered to occur at this temporal position, and the sequences updating is then stopped. The temporal position corresponding to the last element in both side of the sequences define the boundary of the effective part.

When $P=t_0$

$$t_a=t_k^j, ({\sigma }_u>3{\sigma }_d\cap {\epsilon }_{m-k}=\min (\{t_m^u-t_k^d>0\}));$$

(6)

When $P=t_{\mathrm{end}}$

$$t_b=t_k^j, ({\sigma }_u>3{\sigma }_d\cap {\epsilon }_{m-k}=\min (\{t_m^u-t_k^d>0\})).$$

(7)

Based on the previously mentioned selection principle, the effective part $x_2(t_2)$ can be finally determined in the interval $\left[t_a,t_b\right]$.

2.3.2. Noise Smoothing by Preserving Signal Moments

The original echo signal contains a reflected pulse and noise. To suppress the influence of random noise, the waveform data must first be denoised [27].

Let $y(n)$ represent the response function of the discrete system:

$$y(n)=f(n)+(e_f(n)+e_w(n)),$$

(8)

where $n$ is the number of elements in a discrete sequence, $e_f\left(n\right)$ represents the model error in the original signal, and $e_w\left(n\right)$ is the random error caused by white noise. When the echo intensity value has sufficient data, the method of reducing the noise influence while preserving the original signal moment can be adopted [12]. The principle of denoising is to ensure that the signal moment is equal before and after filtering. When the three-order signal moment is preserved, the impulse response is determined as follows:

$$h(n)=3\frac{(3K_1{}^2+3K_1-1-5n^2)}{(2K_1-1)(2K_1+1)(2K_1+3)}\hspace{1em}\left|n\right|\le K_1,$$

(9)

where $2K_1$ represents the order of filtering. We conducted many experiments to estimate the order of filtering and present the test data in this paper. Based on the data from the following experiments, noise rejection was most obvious at $K_1=10$, which was an empirical value.

2.3.3. Digital Differential Low-Pass Filter

A low-pass filter is used to remove or weaken the frequency of the other signals while preserving the reflected frequency from the sea bottom. The key to using this method is determining a method of obtaining the cut-off frequency.

The discrete echo signal represented by the Gaussian function is as follows:

$$y_G(n)=A\exp \left(-\frac{{\left(n-n_{\max }\right)}^2}{2b^2}\right),$$

(10)

where $A$ represents the amplitude, $T$ is the sampling period, $t_{\max }$ is the peak moment, $\sigma$ is the standard deviation of the waveform sequence, $n_{\max }=\frac{t_{\max }}{T}$ represents the sampling rate at the peak of the reflected waveform, and $b=\frac{\sigma }{T}$ is the control parameter of the wave width. Equation (11) can be obtained by Fourier transformation of Equation (10):

$$Y_G(v)=b\sqrt{2\pi }\exp \left(-\frac{b^2{T}^2{v}^2}{2}\right),$$

(11)

where $v$ is the frequency and $b$ represents the width of the reflected echo waveform, which is related to the properties of specific reflectors. If the $b$ value of the bottom in the target area can be determined, the approximate range of the echo signal frequency can be obtained using Equation (11). Related studies have shown that 90% of the types of bottom reflections have values of $b>4.88$. The relationship between $v_c$ and $b$ can be derived as $v_c=0.683/b$; thus, $v_c\le 0.14$ [12].

The low-pass filter is designed as a digital window filter used in the noiseless condition:

$$h(n)=-\frac{1}{\pi }\left(\frac{\sin n\pi v_c}{n^2}-\frac{\pi v_c\cos n\pi v_c}{n}\right),$$

(12)

where $n=1,2,\dots K_2$ and $K_2$ is the order of the filter.

2.4. Gaussian Half-Wavelength Progressive Decomposition

2.4.1. Gaussian Component Parameter Initialization

The received echo signal involves two parts: target reflection and signal attenuation. Under the action of the superposition of the adjacent components and propagation delay, it can be found that the position of the peak in the echo signals is shifted to the left. Based on this situation, Henry Wong proposed that the average value of the temporal position of the local maximum and the minimum of the first derivative can be used instead of the position of the zero, which could effectively eliminate the impact of the peak left shift [12].

$${\mu }_{i\_0}=\frac{1}{2}(t_{pi\_2}-t_{pi\_1}), f^{″}(t_{pi\_2})=f^{″}(t_{pi\_1})=0 \mathrm{and} f^{\prime }(t_{pi\_2})<f^{\prime }(t_{pi\_1})$$

(13)

In Equation (14), $t_{pi\_2}$ and $t_{pi\_1}$ are the temporal positions corresponding to the local extreme values of the first derivative $f(t)$.

By using the method proposed in this paper, the full waveform is decomposed into several Gaussian components according to the time sequence of the signal reception, which reduce the problem of the position shift caused by the waveform overlap. Specific details are discussed later in the Results and Discussion section.

The temporal position of the peak is always a non-integer. The corresponding peak value of $A_{i\_0}$ can be obtained by using the interpolation method based on the non-integer value ${\mu }_{i\_0}$:

$$A_{i\_0}=f({\mu }_{i\_0}).$$

(14)

In the process of gradually decomposing the waveform, the first peak value of the previous residual sequence is selected as the Gaussian component that includes the temporal position and the corresponding amplitude at the current iteration. The key problem is identifying the width of the Gaussian component. Therefore, an iterative method is used to identify the standard deviation ${\sigma }_i$ of the Gaussian function. If one width can satisfy the optimization condition, it is chosen as the wave width of the Gaussian component. Specifically, we use the temporal distance $\mu -t(x_{\sigma })$ as the initial value of the wave width to establish a Gaussian function, where ${\mu }_i$ is the temporal position of the peak and $t(x_{\sigma })$ is the temporal position of the left inflection point. If the waveform value of the corresponding temporal position is less than $m$ times $A_i$ ($A_{im}=m\times A_i$), $\mu -t(x_{\sigma })$ is replaced by $\mu -t(A_{im})$. In addition, $m$ usually takes values of 0.1, 0.3 or 0.5. To guarantee the accuracy of the iterative results, a smaller $m$ could improve the iteration scope.

$${\sigma }_{i\_0}=\{\begin{array}{c}\mu -t(x_{\sigma })\hspace{1em}\hspace{1em}\hspace{1em}{x}_{\sigma }>A_m\hfill \\ \mu -t(A_m)\hspace{1em}\hspace{1em}0<x_{\sigma }\le A_m\end{array}$$

(15)

The three initial parameters $\left[A_{i\_0},{\mu }_{i\_0},{\sigma }_{i\_0}\right]$ that fit the Gaussian function can be initially determined as follows:

$$f(t_i)=A_{i\_0}{e}^{-{(t-{\mu }_{i\_0})}^2/2{\sigma }_{i\_0}^2}.$$

(16)

2.4.2. Component Parameter Optimization and the Termination Condition of the Iteration

The standard deviation ${\sigma }_{i\_0}$ of the Gaussian function determined by the initial value has errors. To obtain an accurate fitting effect, the method of multiple calculations is used to obtain the result.

$${\sigma }_{i\_n+1}={\sigma }_{i\_n}+\Delta \sigma$$

(17)

The value of $\Delta \sigma$ can influence the accuracy of the fitting. If the step length of the progressive calculation is smaller, the fitting process will be finer.

When $\min [f(t_j)-x(t_j)]\ge 0$, the termination condition of the iteration is

$${\sigma }_{v\_ij}=\min (\frac{{\displaystyle \sum _{j=0}^N{[f(t_j)-x(t_j)]}^2}}{N-1}),t_0\le t_j\le {\mu }_i,$$

(18)

where $x(t_j)$ is the fitting value of the half-wavelength of the Gaussian function at the $t_j$ moment and $N$ is the number of sampling points of the half-wavelength range. To ensure the residuals of the fitted waveform are positive and the fitting of the function is satisfactory, the iteration is stopped, when the residual value of the half-wavelength is greater than zero, and its variance ${\sigma }_{v\_ij}$ reaches the minimum.

2.4.3. Waveform Fitting and Parameter Optimization

Since each Gaussian decomposition addresses the residuals in the previous fitting, the residual waveform can be fitted again to restore the original waveform. Additionally, the nonlinear least squares Levenberg-Marquardt (LM) method is used to optimize the parameters of the Gaussian components to improve the fitting accuracy of the full waveform.

With the observation data $(t,y)$, the Gaussian component function is as follows:

$$f_i(t,q)=A_i\exp \left[\frac{-{(t-{\mu }_i)}^2}{2{({\sigma }_i)}^2}\right], q\sim \left(A_i,{\mu }_i,{\sigma }_i\right).$$

(19)

The iterative formula of the LM algorithm is as follows:

$$\widehat{q}=q^{(0)}+{\left[J^T(t,q^{(0)})J(t,q^{(0)})+\lambda E\right]}^{-1}·J^T(t,q^{(0)})\left[y-f(t,q^{(0)})\right],$$

(20)

$$\Vert \widehat{q}-q^{(0)}\Vert <\epsilon$$

(21)

where $J(t,q)$ is the Jacobian matrix of the Gaussian component, $\lambda$ is the damping coefficient, and $\Vert \widehat{q}-q^{(0)}\Vert$ is used as the criterion to judge the termination of iteration. If the value of $\Vert \widehat{q}-q^{(0)}\Vert$ is larger than the threshold $\epsilon$, $\widehat{q}$ would be calculated as the new $q^{(0)}$ until the termination condition (21) is established.

2.5. Selection of the Reflected Waveform

The reflected waveform is not a standard Gaussian function curve, and it is usually characterized by a “trailing” feature that decays over time (Figure 4). Because of the influence of laser pulse propagation in background noise and energy loss, the temporal position of the reflected component is usually erroneously estimated. Therefore, the decomposed results must be selected, merged and denoised to more precisely describe the actual measurements of bathymetry.

Considering the attenuation of the echo pulse energy and the overlap of weak waveforms, the position of the reflected echo should be the moment that the echo energy is mutated during the energy attenuation process. Therefore, based on the Gaussian progressive decomposition, the reflection component in the trailing part is extracted by calculating the second derivative of the component peak sequence. Additionally, based on this method, preselected components are further screened by setting the following limiting condition, and then the final results are obtained.

2.5.1. Intensity

To discriminate the reflected signal from the components, the statistical characteristics of the background noise should be used to eliminate the components that would be influenced by the noise effect or the pulse energy attenuation. If it is assumed that $x_1(t)$ and $x_3(t)$ in Figure 3 are not affected by the reflection, the fluctuation of the residual part is mainly influenced by the background noise, which is usually treated as an additive Gaussian white noise.

Setting three times the standard deviation of the background noise as the confidence interval:

$$X_i\in \left[-3{\sigma }_{noise\_i},3{\sigma }_{noise\_i}\right],$$

(22)

where $X_i$ is the amplitude of the peak. According to the Formula (4), we choose the maximum ${\sigma }_d$ of both sides as the ${\sigma }_{noise\_i}$.

When the intensity value of the decomposed component is in the confidence interval, the intensity value is considered to be derived from the background noise and cannot be treated as a reflection. In other words, the intensity value of the real component must lie outside this interval.

2.5.2. Waveform Width

Because of the attenuation of the signal in the transmission process, the reflected waveform received by the ALB system becomes wider and the intensity decreases. Therefore, the reflected signal in the waveform must be larger than the width of the emission waveform. The attenuation of the signal is related to the intensity of the original emission signal, the propagation distance, the type of medium and other factors, and it is difficult to generalize. The limit conditions for the wave width are as follows:

$$\widehat{\sigma }>{\sigma }_{tr},$$

(23)

where $\widehat{\sigma }$ represents the wave width of the current Gaussian components and ${\sigma }_{tr}$ represents the wave width of the emission waveform. An echo Gaussian function less than the width of the transmitted wave is mainly the result of noise fluctuations; therefore, this peak should be ignored to avoid errors from misattribution.

2.5.3. Relative Position of the Component Peaks

The superposition of multiple echoes occurs during the reflection of the pulse signal. The complexity of the echo wave superposition is mainly caused by the interaction between the amplitude of the component and the echo component interval. After the simulation experiment, Lin suggested that if the interval of the adjacent peak is less than one pulse length, then only one obvious peak is located; however, if the interval is less than a half pulse length, then the temporal position of the different waveforms would be difficult to distinguish [28]. The law of wave superposition is presented in

Table 1. Morphological characteristics and location of adjacent echoes.

Time Interval of Adjacent Components	Morphological Shape of Echo Waveform
$(t_{i+1}-t_i)<\frac{1}{2}{W}_i$	The waveform is similar to one Gaussian function, which can be used only to detect the peak position after superposition.
$\frac{1}{2}{W}_i<(t_{i+1}-t_i)\le W_i$	There is obvious superposition of adjacent components, frequently resulting in loss of peak detection.
$(t_{i+1}-t_i)>W_i$	There is no obvious superposition phenomenon; thus, the different positions of the echo signal can be detected.

.

In the table above, $t_i$ represents the temporal position of the Gaussian peak and $W_t$ represents the length of the pulse. Considering the relative position of the adjacent reflected peaks, the restrictive condition for the interval between the different reflected components should be set as follows:

$$(t_{i+1}-t_i)>\frac{1}{2}{W}_i.$$

(24)

2.6. Metrics for the Comparison

To test the effectiveness and stability of the algorithm, we choose the following metrics for evaluation and comparison:

(1)

The root mean squared error (RMSE) between the estimated depth of the GHPD method and the IWD system is:

$$RMSE=\sqrt{\frac{{\displaystyle \sum _i^M(z_i-z_i^{\prime })}}{M}},$$

(25)

where $z_i$ represents the estimated depths, $z_i^{\prime }$ represents the reference depths and $M$ is the number of waveforms.$RMSE$ can be used to measure the accuracy of the detection.

(2)

The success rate is the percentage of successfully processed full-waveforms; it is given by

$$S=\frac{M_S}{M}\times 100\%,$$

(26)

where $M_S$ is the number of successfully processed full-waveforms that could be detected in more than or equal to two reflections.

(3)

R-squared (R²) represents the fitness of the depth that was successfully detected and is given by

$$R^2=1-\frac{{\displaystyle \sum _i^{M_S}\left(z_i-z_i^{\prime }\right)}}{{\displaystyle \sum _i^{M_S}\left(z_i^{\prime }-\overline{z}\right)}}.$$

(27)

We choose the average of the $z_i^{\prime }$ as the $\overline{z}$ in this paper.

3. Results and Discussion

3.1. Selection of the Effective Part of the Signal

The conventional approach to selecting the effective part of the signal sets a signal intensity threshold as the basis for judgment. Although this method is easy to operate and implement, problems with its effects have been observed:

• Guaranteeing the signal integrity is difficult. Because the method uses the threshold to limit the fluctuation range and mainly depends on the selection of this value, and it is easy to cause discontinuities in the effective part of the signal.
• In the process of effective part selection, a unified and single threshold cannot be implemented because the fluctuations caused by noise interference cannot be considered, and erroneous judgments about the effective interval may occur. In particular, the unified threshold cannot satisfy all the conditions of a noisy environment in areas with deeper water or great changes in the water environment.

The Aquarius full-waveform data are used as an example, and the system implements a single-band blue green laser to scan the experimental area at a scanning frequency of 1 GHz. The number of the discrete time signal in the waveform is 288. According to the method presented in Section 2.3.1, every waveform is extracted to obtain the effective part. Figure 5c shows that the effective part of the echo waveform has better integrity and continuity than the threshold method because of its use of the statistical characteristics of the background noise, and the selection interval exactly contains the reflected effects of the water surface, water column, and bottom.

By selecting the effective part of the waveform, the data redundancy of the intensity sequence was objectively reduced and the computing efficiency of the data processing was directly improved.

Table 2. Efficiency comparison of single-waveform processing.

Section	Average Number of Sampling Points	Average Processing Time of Single Waveform (s)	Time Ratio of Waveform Processing (%)
The effective part	86	0.106	4.93
The complete waveform	288	2.155

shows the statistics of 283,110 full-waveform data points in the experimental area and compares the statistical results before and after the effective part of the waveform was selected. In the process of calculation, the step of the wave width used for the iterative calculation was 0.2. Because the complexity of the waveform is closely related to the bottom topography and the water environment, differences in the number of Gaussian components that could be decomposed from different waveform data were observed, and the data in the table reflect only the average condition of the processing efficiency in the experimental area. The selection of the effective part of a waveform has positive significance for improving the efficiency of batch processing and guaranteeing the quality of the results. Moreover, the selected effective waveform increases the centrality of the target section and avoids the influence of noise and non-reflective signals outside the effective part on the results.

3.2. Results of Waveform Fitting

To test the effect of the method on waveform fitting, a subset of the full-waveform data from Experimental areas A and B was selected as the test object. After smoothing the noise and filtering in the frequency domain, the full-waveform data of the two different bottom materials were decomposed by the GHPD and the parameters of each Gaussian component were obtained. The waveform was then fitted by the superposition of all components. The following figure directly shows the difference between the fitted waveform and the original discrete sampling echo intensity.

Figure 6 shows the fitting effect of the GHPD method intuitively. After the reconstruction of the Gaussian component, the effective part of the waveform is reconstructed. The average of the differences between the original discrete waveform and the reconstructed waveform is 0.1007, the standard deviation is 13.7066, and the Grey Euclid Relation Grade [29] is 0.9214. The fitting waveform after Gaussian component reconstruction displays good fitting accuracy in the effective section of the waveform, and the reconstruction results provide a continuous time function that is consistent with the characteristics of the original discrete waveform data. To further understand the overall effect of regional waveform fitting, we selected portion of the original waveform in Experimental areas A and B for analysis; in each area, the number of selected data points was 7000 (Figure 7).

The depth range in the area was approximately 6–20 m. In this paper, the method of relational analysis was used to evaluate the GHPD method. Because the Grey Euclid Relation Grade overcomes the defect of ignoring the local difference in the calculation, we selected this parameter as the evaluation indicator for the waveform fitting.

The statistical analysis of the relation grade of the fitting result and original waveform showed that the frequency of the relation grade was usually concentrated in the interval $[0.8,1]$ and the relation grade was much greater than the average of the experimental data. This result showed in Figure 7. that the continuous waveform fitted by the GHPD method had high similarity to the original discrete waveform, mainly because the method of progressive fitting is based on the order of echo signal receiving time, which makes the estimation of the initial value of the Gaussian component more reasonable and accurate. Using the LM method to optimize the parameter of each Gaussian component by iterative computation increases the consistency of the superposition of the waveforms of all components with the original waveform, and use of the optimized parameters provides the best fitting effect. Additionally, large amounts of full-waveform data, including much data from complex environmental regions, were used in this experiment for processing and analysis, the results demonstrated that the GHPD method has better algorithm stability in complex propagation environments and under bottom conditions.

3.3. Decomposition of Echo Waveforms and Screening of Reflected Components

The GHPD method uses the same time sequence as the signal reception for decomposing the waveform data, thereby guaranteeing the accuracy of initial value estimation during processing. The reflected waveform was progressively decomposed using the method proposed in Section 2.2 until the entire waveform was processed (Figure 8).

The LM method was then used to optimize the parameters. Several Gaussian components and their parameters are shown in

Table 3. Gaussian component parameters.

Component	Amplitude (${\mathit{A}}_{\mathit{i}}$)	Temporal Position (${\mathit{\mu }}_{\mathit{i}}$)	Gauss Width (${\mathit{\sigma }}_{\mathit{i}}$)
1	109.104	25.159	6.834
2	60.542	33.483	6.531
3	29.684	39.314	3.856
4	33.378	44.950	5.597
5	16.123	50.030	4.120
6	18.467	56.084	5.137
7	283.031	71.697	6.586
8	93.075	77.915	4.658
9	26.639	82.346	3.786

.

Due to energy attenuation and scattering during propagation, the Gaussian components shown in the above table are not all caused by reflection. To obtain the real reflected component, screenings must be performed using the method described above. The specific screening conditions are shown in

Table 4. Limiting conditions for screening the reflected components.

Conditions	Amplitude	Wave Width (ns)	Component Peak Interval (ns)
threshold	>6.10	>4.88	>8

.

Figure 9 shows the results after waveform processing, parameter optimization and reflected component screening. The components shown in red are the reflected components; they show the specific temporal position, intensity and width of the reflected wave. The underwater topography and bottom material were studied by analyzing the reflected parameters.

When two neighboring echo waveforms are close enough together, the overlap between them will occur. The “trailing” of the previous component will affect the peak position of the next waveform component and make it offset (Figure 10). For the waveform processing methods that need to determine the number of wave components and the position of the main peak value first, the overlap will cause a deviation in the temporal position in the depth estimation.

The GHPD method is used to decompose the waveform according to the time sequence of signal reception and selects the reflected components, which can effectively reduce the effect of the peak-position offset. In this section, the echo waveform signal with a depth of 3 m was produced by using the Water LiDAR (Wa-LiD), a waveform simulation tool for LiDAR bathymetry [19,30,31], that the incident angle was 15° and the relative refractive index was 1.333. The simulation waveform isshowed in Figure 11.

The GHPD method and conventional Gaussian decomposition (CGD) method were both used to analysis the waveform data above. Figure 12 directly reflects the difference between the two methods.

In Figure 12d–f, after removing the reflected component of water surface, the bottom parts with great amplitude were selected to decompose using the CGD method given that the bottom parts are still affected by the superposition with the previous scattering waves. Therefore, the characteristic of the bottom reflection is difficult to estimate accurately by the CGD. Through decomposing the waveform according to the time sequence, it can be seen in Figure 12a–c that the influence of the overlapping caused by the previous components has been eliminated using the GHPD, and the reflection components obtained by the GHPD method are more consistent with the original data in Figure 11 than the CGD.

Table 5. Effect comparison of the GHPD and the CGD.

Comparison Index		Simulate Waveform	GHPD	CGD
Surface reflection component	$A$	97.37	101.02	103.38
	$\mu$	49.323	49.419	50.2
	$\sigma$	3.4303	3.4353	3.8219
Bottom reflection component	$A$	16.288	17.057	26.361
	$\mu$	76.519	76.73	75.597
	$\sigma$	3.6068	3.5047	3.3973
propagation time (ns)		27.196	27.311	25.397
Propagation Slant distance (m)		3.0582	3.0712	2.8558
Depth (m)		3	3.0128	2.8015

reflects that the GHPD method could more accurately determine the temporal position of the reflected component in the case of wave overlapping. Meanwhile, the processing of the decomposition is more reasonable, and the estimation of the reflected waveform parameters is closer to the actual situation by using the GHPD method.

3.4. Regional Processing Effect

The seabed point cloud directly reflects the features of the topography. To macroscopically analyze the main features of the GHPD method, we selected the full-waveform data of typical topographies in the experimental area for processing, and the results of regional point-data cloud processing were compared with those of the Aquarius software.

After removing the surface laser points, the shape of the underwater reefs can be clearly seen in Figure 13b,d, whereas the laser foot points of the underwater terrain in Figure 13a,c were buried in the cloud of multiple echoes and were difficult to separate due to fluctuations in the water surface and reverberations in the water. In the experiment, IWD detected more reflection points but was also more susceptible to surface fluctuations, water impurities, etc. In the experimental area, the GHPD method was applicable to shallow areas in which the water depth was less than 5 m. Although the number of reflected echoes that could be determined by the GHPD method was less than the number determined using the system software, the submarine topography features were all intact.

Because of the screening mechanism, the GHPD method has a certain inhibitory effect on the propagation noise caused by water surface fluctuations and water impurities.

3.5. Consistency with the Data Processing Software

To verify the effect of the GHPD method on underwater terrain, the partial data at depths of approximately 2 m~16 m in Experimental area A and Experimental area B were selected for analysis in this section. The seabed laser point cloud could be obtained by the full-waveform processing with refractive index correction. Figure 12 shows the seabed laser point cloud derived using the two methods in the different experimental areas.

By comparing the processing results of the two methods in different experimental areas, it can be intuitively found that the GHPD algorithm and IWD show little difference in the seabed morphology, and the coverage of the point cloud is relatively similar, which verifies the batch processing results of this method are consistent with that of mature commercial software.

It can be seen from the local magnified part of Figure 14 and

Table 6. The overall number and the local density of laser points.

Experimental Area	Area (${\mathbf{m}}^2$)	Method	Number of Laser Points	Density (/${\mathbf{m}}^2$)
A	91103.02	GHPD	144936	1.59
		IWD	165230	1.81
B	196716.82	GHPD	309411	1.57
		IWD	342657	1.74

that the density of the point cloud retrieved by using GHPD method is slightly inferior to that of the IWD system, especially in the deeper areas or the scanning edge, but the density of the results achieved using the two methods could be satisfied to the requirements for the general bathymetry. The GHPD method can extract more information, like the width of the reflection waveform, which could be used to determine the feature of the reflected surface or the backscatter.

To analyze the precision and stability of the algorithm, the partial full-waveform data are selected and processed, sorted and compared according to different depth intervals. The results processed by IWD, which is the stable commercial software, are taken as reference values. The relevant statistical metrics were then calculated as shown in the

Table 7. Performance of the GHPD method in Experimental area A.

Depth Interval (m)	Number of Detectable Waveforms	RMSE (m)	S (%)	R2
2~4	15,306	0.052	83.122	0.904
4~6	32,191	0.059	80.299	0.991
6~8	37,630	0.056	93.141	0.992
8~10	35,710	0.053	92.358	0.992
10~12	20,046	0.050	89.389	0.992
12~14	4040	0.045	87.946	0.988
14~16	13	0.033	88.769	0.688
Mean	Total: 144,936	0.050	87.718	0.935

and

Table 8. Performance of the GHPD method in Experimental area B.

Depth Interval (m)	Number of Detectable Waveforms	RMSE (m)	S (%)	R2
2~4	1551	0.050	82.721	0.972
4~6	7949	0.050	90.937	0.979
6~8	88,805	0.043	88.433	0.991
8~10	124,393	0.048	90.364	0.992
10~12	60,869	0.045	94.084	0.993
12~14	25,834	0.041	95.545	0.992
14~16	10	0.060	90.000	0.977
Mean	Total: 309,411	0.048	90.298	0.985

and Figure 15.

The waveform processing results show that the processing accuracy, the method success rate and the stability of the GHPD method in Experimental area B (coral reefs and coral reef debris) are on the whole slightly better than those in Experimental area A (fine sand). The waveform data in this section mainly covers the depth range from 2 m to 16 m. Except in the 14–16 m range, where the water is too deeper and the effective echo that the Aquarius system can obtain is less, there is no significant change in the comparative metrics at other depth intervals. The result shows that the GHPD method has good stability in estimating the certain depth of water. The average success rate of each depth interval waveform is greater than 87%, and the RMSEs of Experimental areas A and B are not more than 0.06 m; thus, the results fully satisfy the specifications of the International Hydrographic Organization (IHO) for bathymetry [32]. This result proves that the results obtained using this method can be applied to actual bathymetry. The success rate of the GHPD method is lower in areas of this type where the experimental depth is too shallow (<2 m) or too deep (>14 m). This problem indicates that the GHPD method is worth further considering in terms of the signal extraction of such feature regions. The data in Figure 15 show that the metrics are better when the depth is approximately 10 m, indicating that the accuracy of detection of the GHPD method is optimal in this depth range.

4. Conclusions

To resolve the problems associated with full-waveform airborne laser bathymetry data processing, a full-waveform data processing strategy based on GHPD was proposed, and the processing flow of this method was expounded in detail. The results obtained with the GHPD method were consistent with those obtained with commercial software. The effectiveness of the GHPD method in full-waveform processing and information extraction was proven, and the following conclusions were reached:

• In the data preparation stage, the effective part of the full-waveform data must be selected and smoothed. By using the statistical characteristics of the background noise in the signal propagation process, the effective part of the full-waveform data was extracted as the main object of the study. The experimental results showed that the proposed method can avoid unnecessary noise interference and improve the computational efficiency. Therefore, smooth processing of waveforms was achieved using the method of preserving signal moments to filter the original waveform and considering the reflection properties of the water environment and the statistical characteristics of the waveform data set.
• Based on the waveform decomposition, the GHPD method can simultaneously realize the accurate fitting of full-waveform data. The advantage of the GHPD algorithm is that the fitting and decomposition of the waveform can be realized by Gaussian functions in the time sequence according to the morphological characteristics of the original waveform.
• A proper selection mechanism must be employed in the GHPD method to obtain the reflected component. The screening mechanism can be established on the basis of an analysis of the changes in the roles of signals in the reflection process and the characteristics of pulse intensity, wave width and time interval. After decomposition, all the Gaussian components were judged and selected, and then the corresponding component of the target reflection was obtained. This method objectively reduces the attenuation of the signal in the propagation process and interference of the echo under the influence of the background noise to provide accurate temporal position information that can be used to obtain the target distance.

Additionally, the results obtained using the GHPD method and system-matching software under different substrate conditions were compared by performing a data processing experiment involving a large number of regional full-waveform airborne laser bathymetry data points. The GHPD method and the intelligent waveform digitizer (IWD) algorithm of mature commercial software presented consistent estimates of the laser propagation distance in water. Additionally, the two methods were similar with respect to the water depth detection, and the difference in the results between the two methods was small.

Nevertheless, there are still many shortcomings in this calculation strategy, such as the calculation complexity, the selected correlation threshold and parameters, empirical judgment, and the problem of self-adaptation. However, the proposed method can still improve the performance of airborne full-waveform data processing and is worthy of further exploration and research.