Next Article in Journal
Simplifying Sample Preparation for Soil Fertility Analysis by X-ray Fluorescence Spectrometry
Next Article in Special Issue
Remote Sensing Applications in Coastal Areas
Previous Article in Journal
A Gyroscope Signal Denoising Method Based on Empirical Mode Decomposition and Signal Reconstruction
Previous Article in Special Issue
Normalized Method for Land Surface Temperature Monitoring on Coastal Reclaimed Areas
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

A Depth-Adaptive Waveform Decomposition Method for Airborne LiDAR Bathymetry

1
The Institute of Geospatial Information, Strategic Support Force Information Engineering University, 62 Science Road, Zhengzhou 450001, China
2
Science and Technology on Near-surface Detection Laboratory, 160 Tonghui Road, Wuxi 214035, China
*
Author to whom correspondence should be addressed.
Sensors 2019, 19(23), 5065; https://doi.org/10.3390/s19235065
Submission received: 14 October 2019 / Revised: 14 November 2019 / Accepted: 15 November 2019 / Published: 20 November 2019
(This article belongs to the Special Issue Remote Sensing Applications in Coastal Areas)

Abstract

:
Airborne LiDAR bathymetry (ALB) has shown great potential in shallow water and coastal mapping. However, due to the variability of the waveforms, it is hard to detect the signals from the received waveforms with a single algorithm. This study proposed a depth-adaptive waveform decomposition method to fit the waveforms of different depths with different models. In the proposed method, waveforms are divided into two categories based on the water depth, labeled as “shallow water (SW)” and “deep water (DW)”. An empirical waveform model (EW) based on the calibration waveform is constructed for SW waveform decomposition which is more suitable than classical models, and an exponential function with second-order polynomial model (EFSP) is proposed for DW waveform decomposition which performs better than the quadrilateral model. In solving the model’s parameters, a trust region algorithm is introduced to improve the probability of convergence. The proposed method is tested on two field datasets and two simulated datasets to assess the accuracy of the water surface detected in the shallow water and water bottom detected in the deep water. The experimental results show that, compared with the traditional methods, the proposed method performs best, with a high signal detection rate (99.11% in shallow water and 74.64% in deep water), low RMSE (0.09 m for water surface and 0.11 m for water bottom) and wide bathymetric range (0.22 m to 40.49 m).

1. Introduction

Airborne LiDAR bathymetry (ALB) is a technique for measuring the depths of moderately clear, near-shore coastal waters and lakes with a high-powered, high-pulsed green laser from a low-altitude aircraft (200–500 m above ground level (AGL)) [1,2,3]. This technique can provide high-density, high-accuracy, and three-dimensional bathymetric data safely and with high efficiency compared with shipborne sonar. The swath width of an ALB system roughly ranges from 0.5–0.75 times AGL, namely 100–400 m, indicating that the area covered within 1 h ranges from about 20–60 km2 [4]. Since ALB is cost-effective and time-saving, it has been widely used in coastal mapping, sediment budgets and seabed protection [5,6].
Since the transmission of green laser involves air and water, the waveform processing for the ALB system sometimes could be challenging. Water absorption and scattering influence the extraction of the water surface and bottom returns from the waveform of a green channel, especially for turbid water and shallow water (a depth lower than 0.5 m) [7,8]. Water column scattering, which is much more powerful than atmospheric scattering, leads to a low signal-to-noise ratio (SNR) in the received waveform. Two main problems in the full waveform processing for ALB are given below:
  • Mixed peaks in surface return: The waveforms are the convolution of the emitted pulse and the target cross section, and are digitized by the receiver. The limited full width at half maximum (FWHM) of the emitted pulse and sampling rate in the LiDAR digitizer induce the peak stretching, leading to a mixed peak of surface return and water column scattering. Especially when the water is extremely shallow, the bottom return will also be included in the mixed peak. Taking the mixed peak as the surface return may introduce errors ranging from 10 cm to 25 cm [9].
  • Weak bottom return in deep or turbid water: The pulse energy decreases exponentially with depth in the water column, and the decrease rate is positively associated with water turbidity, resulting in a rather weak bottom return in deep or turbid water [10].
For the first existing problem, some ALB systems use an NIR channel or a Raman channel to measure the surface return. Pe’eri and Philpot [11] studied the relationship between the shape of the Raman waveform and the water depth and used the Raman channel to measure the water depths (shallower than 2 m). However, this algorithm may vary with water clarity. Allouis et al. [12] found that the Raman signal is not reliable for depth estimation, since the Raman signal is sensitive to water characteristics and proposed a depth estimation method using NIR and green fitted waveforms. Another alternative is to establish a water surface model to determine the position deviation between the mixed peak and the water surface, so that the single green laser ALB systems can be used to estimate the water depth. Mandlburger et al. [9] analyzed the near water surface penetration (NWSP) properties of green laser signals. Zhao et al. [13] studied the factors that influence the NWSP of a green laser and proposed an NWSP modeling method. This method needs an auxiliary NIR laser and some SSC sampling stations, which sometimes are difficult to obtain, and so waveform decomposition may be the optimal solution to this problem. For topographic LiDAR, the Gaussian model is sufficient for most applications [14]. However, the Gaussian decomposition method is not suitable for ALB because the water column component in the waveform cannot be easily fitted by Gaussian functions [15,16]. The triangular function [17] and quadrilateral function [18] were both introduced to water column fitting but were only verified by simulated data. Ding et al. [19] proposed an improved quadrilateral model for the water column fitting which shows a better fit to the field data compared with the quadrilateral function. For very shallow water, a surface-volume-bottom (SVB) algorithm was proposed by Schwarz et al. [20] and was applied to measure a riverbed. However, waveforms at different depths vary greatly, and a single model is only suitable for waveforms in a specific depth range.
For the second problem, the key is detecting the bottom signal from low-SNR waveforms. Numerous waveform processing methods have been developed for waveform de-noising or signal enhancement. Saylam et al. [21] used a moving-average filtering algorithm for waveform smoothing. Pan et al. [7] introduced a continuous wavelet transformation (CWT) to project the signal into a continuous time and scale subspace and found that more signals can be detected from the reconstructed waveform while many signals are still undetectable. Wang et al. [15] examined some waveform processing methods and concluded that the Richardson–Lucy deconvolution (RLD) [22] is good at dealing with waveforms with a very shallow depth and a weak bottom response, and the average square difference function (ASDF) [23] can better cope with noise. Richter et al. [24] proposed an attenuation correction procedure to improve the detectability of water bottom signals. Launeau et al. [25] proposed a waveform processing method including smoothing and edge enhancing to increase the detection rate of the bottom signal in moderately turbid water. The de-noising and signal-enhancing methods can improve the detection rate, but the detection accuracy is still limited by the waveform sampling interval, which may also require waveform decomposition.
Although the previous methods have partially solved the above problems, there is no single best waveform processing strategy across all applications [7,26]. The received waveform influenced by the water depth and the optical attenuation of the water column is varied and complicated. If each component of the waveform is reconstructed, then both the water surface and bottom signals can be precisely detected, and the accurate water depth can be estimated. In this study, we proposed a depth-adaptive waveform decomposition method for green channel waveforms of ALB. Waveforms are first classified by depth based on the shape. In waveform decomposition, two models are developed for waveforms with different depths based on the calibration waveform and an improved quadrilateral model, and a trust region algorithm is introduced to solve the model parameters. Field data from an airborne topo-bathymetric LiDAR system and simulated data are used to assess the accuracy of the detected water surface and bottom positions. The rest of the paper is organized as follows. The field datasets used in this paper are introduced in Section 2, followed by the detailed description of the proposed wave decomposition method. The experimental results are presented and discussed in Section 3. Finally, Section 4 concludes this paper.

2. Materials and Methods

2.1. The Mapper5000 System and Study Area

The field data were collected by an airborne topo-bathymetric LiDAR (Mapper5000) system in the Qilianyu Islands, Hainan Province, China. This ALB system was developed by the Shanghai Institute of Optics and Fine Mechanics, Chinese Academy of Science. The system performed an elliptical scan at a scan angle of ±15° giving a field of view (FOV) of 30° over a range of 16.8 km × 2.2 km, and the detailed parameters are shown in Table 1. The sensor consists of four channels, including one NIR channel (1064 nm, avalanche photodiode (APD)) and three green channels (532 nm, photomultiplier tube (PMT)). These three green channels are set differently. PMT1 and PMT2 have a shallow field angle, while PMT3 has a wide field angle. The instantaneous field of view (IFOV) of NIR, PMT1 and PMT2 is 6 mrad, and the IFOV of PMT3 is 6‒40 mrad. The receiving directions of PMT1 and PMT2 are perpendicular to each other. Because signals in PMT1 waveforms are generally stable, this study focuses on the processing of PMT1 waveforms, and waveforms in other channels are used as references. Besides this, the calibration waveforms collected in the laboratory are utilized as the approximation of the emitted signal. With the detected bathymetric signals, point clouds can be generated by a data processing software. This data processing software is designed for the Mapper5000 system including geo-calibration and refraction correction for the point clouds. Detailed descriptions of the Mapper5000 system and its point clouds generation software are appended in Appendix A and Appendix B, respectively.
For the statistical analysis of waveforms, we choose four bathymetric points (P1–P4) with different depths. In order to access the performance of the proposed waveform decomposition method, especially in areas where the water depth is deep or shallow, two representative datasets, Dataset 1 and Dataset 2, were selected from two different strips. Figure 1 shows the locations of the strips, the datasets and bathymetric points in the study area.

2.2. Workflow

The original PMT1 waveforms are processed with the depth-adaptive waveform decomposition method to extract the accurate positions of the bathymetric signals. As shown in Figure 2, this is a multistep process including preprocessing, signal detection and waveform decomposition. The first step was to determine the useful range of the waveforms, classify the waveforms by depth and improve the signal resolution by RLD. After deconvolution, the second step was to detect the signals with an adaptive threshold, and the results were used as the initial values of the following waveform decomposition. PMT1 waveforms were classified into two categories, “shallow water (SW)” and “deep water (DW)”. In the last step, the fitting model was selected based on the waveform category. With the initial values provided by the second step, the model parameters were solved by the TR algorithm. After waveform decomposition, the surface and bottom signals were extracted from the waveforms. To evaluate the accuracy of the results, signals detected from NIR waveforms and PMT3 waveforms were used as reference data.

2.3. Preprocessing

2.3.1. Useful Range

To ensure the ALB system can measure both land and water, the system normally records thousands of samples for the received waveform, while the signals only exist in 0.8–5% of it. Thus, if the range of the signals, which is known as the “useful range”, can be determined, the waveform processing will be more efficient. Jutzi and Stilla [27] proposed a criterion that if the waveform is three times higher than the noise power for at least 5 ns, a signal will be assumed to have been found. In this paper, we choose the last 10% of the waveform to estimate the noise. According to the maximum measurable depth and flight height of the system, signals cannot exist in this range. The truncation noise NT and the noise power NP can be estimated using the minimum amplitude and the standard deviation of the waveform in this range, respectively. The noise level NL can be expressed as follows:
NL = NT + 3NP.
By searching the waveform using a criterion that a waveform higher than NL and lasting for no less than 5 ns may indicate a signal, the useful range can be determined as shown in Figure 3.

2.3.2. Waveform Classification

Received waveforms of the green LiDAR are varied with depths as shown in Figure 4. To fit the waveforms with different depths, waveforms are grouped into two categories, “shallow water” and “deep water”, according to a defined parameter S. Because the depths of the water that results in the surface and bottom signal overlap vary with different ALB systems or measured area, it is difficult to classify waveforms with an absolute water depth threshold. Based on the waveform analysis, we found that the shape of the water column scattering in the received waveform is almost fixed in a small measured area (See Section 3.1.1). If the water column scattering can be captured in the waveform, it also indicates that there is no overlap between the water surface and the bottom signal. Thus, the water column scattering can be used as a sign to classify the waveforms. In this study, shallow water and deep water are not distinguished by absolute water depth but by waveform shape. If the surface and bottom signal overlap in the waveform, it is defined as “SW”; otherwise, the waveform is labeled as “DW”.
The parameter S is determined by the following equations:
S = min { R ( t ) } ,
R ( t ) = 1 N m = 1 N [ w C ( m τ ) w R ( m τ + t ) ] 2 ,
where wR is the received waveform, wC is a section of the water column scattering truncated from the received waveforms, N is the sampling number of wC, and τ is the sampling interval of the system. Because the received waveforms in extremely deep water have complete water column scattering (see Figure 4c), wC can be easily obtained from them. Specific details are given in Section 3.1.1. For the “SW” waveform, the water column scattering is embedded by surface and bottom signals, and so the received waveform is less similar to wC, as shown in Figure 5a. Because water column reflectance varies little in a small measured area, the “DW” waveform will have a small S, as shown in Figure 5b. Thus, received waveforms in which S is less than a threshold (TS) are classified as “DW”, while the remaining waveforms are “SW”.

2.3.3. Richardson–Lucy Deconvolution

Waveforms are processed by RLD to increase the signal resolution. Considering that most of the waveforms are likely to have mixed peaks, especially for the first peak which may contains surface and column return, this overlap mostly results from a convolution of the emitted waveform and the target cross section. Thus, a deconvolution process is needed before detecting the signals.
RLD is an iterative deconvolution method which was developed by Richardson [28] and Lucy [29] to recover a blurred image with a known point spread function (PSF). Many studies [22,30,31] have applied this to LiDAR waveform processing in the time domain to estimate the target cross section, where the emitted waveform is regarded as a PSF. The ith iteration can be calculated as
p ¯ i + 1 ( t ) = p ¯ i ( t ) ( w T ( t ) w R ( t ) ( w T ( t ) p ¯ i ( t ) ) ) ,
where “∗” is the convolution product, wT(t) is the emitted waveform which is approximated by the calibration waveform in this paper, wR(t) is the received waveform, and p ¯ ( t ) i + 1 is the ith iteration solution of the target cross section. The initial value of p(t) for the iteration is simply set as wR(t), which has negligible effect on the result. Determining the end point of the iteration is the key of RLD because the noise amplification is due to overfitting. The stopping criterion of RLD used here is detailed discussed in [32].
Figure 6 shows that RLD can improve the signal resolution whether in shallow water or deep water waveforms. In shallow water waveforms, RLD can shorten the FWHM of the surface and bottom signals (see Figure 6a). In deep water waveforms, RLD can partially remove the background noise while keeping the weak bottom signal (see Figure 6b).

2.4. Signal Detection

A number of signal detection methods such as threshold, center of gravity and maximum are often used in conventional LiDAR waveform processing [33]. However, most of them may not work for wR(t) in ALB. We tested some of the conventional signal detection methods, as shown in Figure 7. It could be determined that the center of gravity method may not be suitable for ALB signal detection, due to the component of water column scattering in the received waveform. The threshold and maximum methods can give more reliable results with an appropriate threshold. However, the appropriate threshold may vary in different waveforms.
In this section, an adaptive threshold was used for ALB signal detection instead of a fixed one, as shown in Figure 8. The adaptive threshold T(t) is defined as
T ( t ) = { max { w C } + 3 × N P ,   t < t S w C ( t t S ) + 3 × N P ,   t t S ,
where tS is the minimum point corresponding to the minimum value S of R(t). As the water depth changes, T(t) will be adjusted reasonably, effectively avoiding the influence of high-intensity fake signals on signal detection and improving the reliability of the results.
Candidate signals in deconvoluted waveforms are detected with the maximum method and filtered by T(t). If there are more than two remaining signals, the two signals with the largest amplitude are retained.

2.5. Waveform Decomposition

There are three basic steps in waveform decomposition, including modeling, initialization and fitting. The waveform processed here is denoted by w(t), which is the original received waveform truncated at noise level NL.

2.5.1. Modeling

The key to waveform decomposition is to build a reasonable model. Some models for the ALB waveform have been introduced, such as two Gaussian functions [12], a combination of a Gaussian function, a triangle function and a Weibull function [17], a combination of two Gaussian functions and a quadrilateral function [18], and a chain of exponential segments [34]. Considering that water depth is an important factor influencing the shape of waveform (as shown in Figure 4), two models are proposed for “SW” and “DW” waveforms, respectively.
The fitting model fW(t) can be expressed as
f W ( t ) = f S ( t ) + f B ( t ) + f C ( t ) .
For the surface return model fS(t) and bottom return model fB(t), the existing models use a specific function, such as Gaussian, exponential and harmonic function, to model the surface and bottom return, but are not exact in practice because of the distortion caused by the sensor response, as shown in Figure 9. Thus, we do not use a model with some defined functions, but directly use a transformation C of the calibration waveform to better fit the surface and bottom return:
C ( A , μ , σ , φ ) = A φ ( t μ σ ) ,
where A is the amplitude scaling factor, μ is the time shift factor, σ is the time scaling factor, and φ is the normalized calibration waveform. The calibration waveform used here can be collected in the laboratory or substituted by a waveform received from bare ground, and φ is estimated using the smoothing spline method, normalized in amplitude, and shifted in time to locate the peak at the zero point.
The surface return model fS(t) is given by
f S ( t ) = C ( A S , μ S , σ S , φ ) ,
where AS, μS and σS are the amplitude scaling factor, time shift factor and time scaling factor in function C, respectively.
The bottom return model fB(t) is given by
f B ( t ) = C ( A B , μ B , σ B , φ ) ,
where AB, μB and σB are the amplitude scaling factor, time shift factor and time scaling factor in function C, respectively.
For “SW” waveforms, the water column scattering model fC1(t) is equal to
f C 1 ( t ) = C ( A C , μ C , σ C , φ ) .
where AC, μC and σC are the amplitude scaling factor, time shift factor and time scaling factor in function C, respectively. For “SW” waveforms, the proposed model fW(t) is based on the calibration waveform and is therefore named the empirical waveform model (EW).
For “DW” waveforms, fC(t) is defined as
f C 2 ( t ) = { exp ( f b 2 + g b + h ) ( t a b a ) a < t b exp ( f t 2 + g t + h ) b < t c exp ( f c 2 + g c + h ) ( d t d c ) c < t d 0 e l s e ,
where a, b, c and d are the horizontal coordinates of four boundary points in fC1(t), as shown in Figure 10, and f, g and h are coefficients related to water column scattering. Here, an exponential function with a second-order polynomial is proposed to improve the quadrilateral model presented in [19]. Hence, this model is named the exponential function with second-order polynomial model (EFSP).
Thus, the proposed model fW(t) can be denoted by fW(t, γ) with the unknown parameter vector
γ = ( A S , A B , A C , μ S , μ B , μ C , σ S , σ B , σ C ) ,
or
γ = ( A S , A B , μ S , μ B , σ S , σ B , a , b , c , d , f , g , h ) .

2.5.2. Initialization

Initialization is a significant process in waveform decomposition. Waveform decomposition is a nonlinear nonnegative least-squares problem, which can be solved by many algorithms, but global convergence cannot be promised, especially when there are a large number of parameters to solve.
The first step of initialization is to calculate the initial value γ0. The results of signal detection in Section 2.4, namely the rough positions of the surface and bottom signals denoted by tS0 and tB0, are used for calculating the initial values of the parameters in γ:
( A S 0 , A B 0 , A C 0 , μ S 0 , μ B 0 , μ C 0 , σ S 0 , σ B 0 , σ C 0 ) = ( w ( t S 0 ) , w ( t B 0 ) , 0.5 × w ( t B 0 ) , t S 0 , t B 0 , 0.5 × ( t S 0 + t B 0 ) , 1 , 1 , 1 ) ( a 0 , b 0 , c 0 , d 0 ) = ( t S 0 0.5 × t L , t S 0 + 0.5 × t R , t B 0 0.5 × t L , t B 0 + 0.5 × t R ) ,
where tL is the length of φ to the left of the peak, and tR is the length of φ to the right of the peak. If only surface signal is detected, the value of tB0 will be set as tS0 + 1/2 × tL.
In the EFSP model, the last three parameters in γ are estimated by a simple linear fitting:
A fit = [ t 4 t 3 t 2 t 3 t 2 t t 2 t 1 ] b fit = [ ln ( w ( t ) ) t 2 ln ( w ( t ) ) t ln ( w ( t ) ) ] [ f 0 g 0 h 0 ] T = A fit 1 b fit ,
where the fitting range is set at [tS0 + tR, tB0tL], because the surface and bottom returns have a negligible effect on this part of the waveform.

2.5.3. Fitting

Regarded as a non-linear least-squares problem, the step of fitting is dealt with by many mature theories and algorithms, such as the Gauss–Newton algorithm (GN), Expectation-Maximization algorithm (EM) [35], Levenberg–Marquardt algorithm (LM) [36,37] and Reversible Jump Monte Carlo Markov Chain (RJMCMC) [38]. However, GN, the most traditional method, easily converges at a local optimal solution, and its modification, LM, although it has a global convergence to some degree, may be influenced by the initial value in practice. Thus, the TR algorithm is introduced in this study. Unlike the other algorithms mentioned above, TR can solve a problem with constraints, indicating that the parameters will be calculated in a reasonable range. Furthermore, TR is not a line search algorithm which is carried out along a search direction in each iteration; instead, TR searches the next iterated point in a trust region, which is a neighborhood of the current iterate point. The principle is described below [39].
The cost function Q(γ) is expressed as:
min   Q ( γ ) = i = 1 n [ w ( t i ) f W ( t i , γ ) ] 2 ,   γ R m ,
where n is the number of sampling points in the useful range, and m is the number of parameters in γ. At the kth iteration, Q(γ) is expanded in Tailor at the current iteration point γ (k) with the second-order terms preserved
Q ( γ ) Q ( γ ( k ) ) + Q ( γ ( k ) ) T ( γ γ ( k ) ) + 1 2 ( γ γ ( k ) ) T 2 Q ( γ ( k ) ) ( γ γ ( k ) ) .
We plug d = γγ (k) into Equation (17) and get the quadratic form:
φ k ( d ) = Q ( γ ( k ) ) + Q ( γ ( k ) ) T d + 1 2 d T 2 Q ( γ ( k ) ) d .
The trust region in the current iteration can be expressed as ‖d‖ ≤ rk, where rk is the trust region radius. As the range of γ is given, we can limit the solution in a reasonable range by setting rk. Thus, the cost function in Equation (16) will be translated into solving the following trust region subproblem:
min   φ k ( d ) s . t .   d r k .
By solving Equation (19) using a line searching algorithm, the optimal solution d(k) is obtained. The key part of a TR algorithm is how to judge d(k); that is, whether to accept the current improvement d(k) and how to change the trust domain radius rk. According to the chosen strategy, the TR method has many forms [40]. This paper only uses one of them; that is, the correctness of d(k) is judged according to the ratio of the actual decrease of the function value to the predicted decrease [41]:
ρ k = Q ( γ ( k ) ) Q ( γ ( k ) + d ( k ) ) Q ( γ ( k ) ) φ ( d ( k ) ) .
If ρk ≤ 0, then d(k) approximation fails, and we let γ (k+1) = γ (k). Conversely, if d(k) approximation succeeds, we let γ (k+1) = γ (k) + d(k) and calculate the new trust domain radius rk+1 according to ρk. The new starting point γ (k+1) and the trust region radius rk+1 are re-substituted into Equations (18) and (19). We repeat these steps until the result converges.

3. Results and Discussion

3.1. Experiment Ⅰ: Waveform Classification

3.1.1. Statistical Analysis of Waveforms with Different Depths

Water depth is one of the most important influence factors of the waveform. To classify the waveform by depth reasonably, waveforms of four bathymetric points (P1–P4) with depths around 2 m, 10 m, 20 m and 30 m are randomly selected for analysis (100 waveforms for each bathymetric point). For waveform registration, each waveform was shifted in time to locate the peak of the surface signal at the zero point. The mean and variance curves of the four bathymetric points were calculated as shown in Figure 11.
As depicted in Figure 11a, the bottom signal is most affected by the water depth, and its amplitude shows an obvious decrease with the increase in water depth. Some differences can also be found in the surface signals due to the change of the measurement time and position. However, the amplitudes of water column scattering are almost constant at different depths. The variance curves in Figure 11b show the difference in waveforms at the same water depth. It can be found that the surface and bottom signals are still changing even at a constant water depth. In contrast, the change in water column scattering is negligible.
When the water depth is 30 m and the time offset is 10 ns, the variance quickly drops to 60, indicating that the waveform after this point is slightly affected by the surface signal. When the time offset is equal to 30 ns, the amplitude of the mean curves is below 5% of the intensity of the surface signal; this suggests that the shape of water column scattering in the received waveform is almost fixed in a small measured area, which can be used as a sign to classify the waveforms. Therefore, the mean curve with a time offset between 10 ns and 30 ns can be truncated as wC (see Section 2.3.2).

3.1.2. Distribution of S.

According to the width of the emitted signal, the received waveforms can be divided into three cases:
  • Case 1: Water column scattering is completely covered by the surface and bottom signals.
  • Case 2: The length of water column scattering that is not covered is less than the length of wC.
  • Case 3: The length of water column scattering that is not covered is greater than the length of wC.
In this experiment, waveforms in different cases were selected and S is calculated according to Equations (2) and (3). As shown in Figure 12, the S values of the waveforms belonging to Case 1 are above 4500, while the S values in Case 3 are below 3500. Thus, the threshold TS can be set between 3500 and 4500. In this paper, TS is set to 4000.

3.2. Experiment Ⅱ: Performance Analysis for the Processing Algorithms in Shallow Water

3.2.1. Reference Data

The bathymetric environment is complex and changeable. Figure 13 shows the scope of Dataset 1. It can be seen that the edge of Dataset 1 is adjacent to the island’s coastline, which is the most affected area by tides. Therefore, in accuracy analysis, the reference data must be acquired simultaneously with the field data. However, the bathymetric sonar is inefficient in shallow waters compared to ALB; it is difficult to ensure these two bathymetric methods are performed simultaneously.
To assess the signal detection accuracy, signals detected from NIR waveforms were taken as references. Because the NIR signal does not penetrate into the water column, it can provide an accurate position of the water surface [9]. For shallow water, its waveform in the PMT1 channel has a strong bottom return. Errors occur only when a mixed peak is detected. Thus, the accuracy of surface signal detection can reflect the accuracy of the water depth.
However, NIR signals are not always reliable (sometimes buried below the noise level due to strong absorption and sometimes saturated, see Figure 14). In this experiment, NIR signals are filtered according to the signal strength before being used as references.

3.2.2. Surface Signal Detection

To access the performance of the proposed algorithms in shallow water, the traditional maximum method was applied to analyze the effect of the adaptive threshold (see Section 2.4). Signals in waveforms processed by RLD are both detected by the maximum method with a fixed threshold and adaptive threshold. The corresponding point clouds are referred to as RLD_M and RLD_A. The fixed threshold is set to NT + 3NP. The Gaussian model proposed in [12] is used to process waveforms in very shallow water. To compare the applicability of the model, waveforms are both fitted with the Gaussian and EW model. Their initial values are calculated according to RLD_A, and the generated point clouds are referred to as GD (Gaussian decomposition) and EW. The point cloud only processed by the tradition maximum method is also presented for comparison (denoted by Max).
The water surface points from all the point clouds are compared in height with the reference data (See Figure 15). Most of the inaccurate surface points are lower than the reference data, because the bottom signals are detected as surface signals by mistake or the peak of the surface signal is shifted backward due to the overlap with the bottom signal. The Max is most affected by this problem. Since RLD can improve the signal resolution, this problem has been significantly improved in RLD_M and RLD_A. EW shows a good consistency with the reference data, and most of the errors are within 0.15 m, indicating that the errors of the detected signal position are within one sampling interval of the waveform.
Figure 16 shows the water surface on the selected profile (see Figure 13) detected by the above algorithms (only when both surface and bottom signals are detected will the surface point be displayed in this figure). The water depths on this profile are shallower than 1 m. The height of the water surface should float at 0 m, just as with the reference data. The heights of most points in Max are between −1 m and −0.5, indicating that mixed peaks were detected. In RLD_M, more points have been detected, and the number of error points is significantly reduced, proving that RLD can effectively improve the signal resolution. Although the number of points in RLD_A is less than in RLD_M, most of the missing points are outliers. There are still many outliers in GD. Compared with the other results, the distribution of the points in EW is most similar to the reference data, and the difference is no more than 0.5 m.
Table 2 provides the statistical parameters of the above results, where Dr represents the detection rate, which is defined as the percentage of surface points with an error less than 0.3 m; RMSE denotes the root mean square error of the surface points; min(d) corresponds to the detected minimum depth; and Std. means the standard deviation of the heights of the surface points.
The detection rate of EW is 99.11%, which is the highest. The detection rate of RLD_M, RLD_A and GD is close, but the result of RLD_A is more reliable, which can be found by comparing RMSE. GD and EW have the minimum min(d), followed by RLD_M and RLD_A. Apart from the RMSE, Std. can also reflect the accuracy of signal detection. Since sea level changes very little in a small area, the smaller the Std. is, the more accurate the surface points are. In general, EW performs best, with the highest detection rate, minimum error and maximum bathymetric range. The reference data has a low detection rate (58.54%) and a small Std (0.0839), which indicates that NIR signals are undetectable sometimes but accurate.

3.2.3. Adaptability Analysis of the EW Model

To further study the effects of different models on waveform decomposition, three representative waveforms are selected from Dataset 1. The corresponding depth of the first waveform is more than 1 m, and the initial values provided by RLD_A are exact. The depths of the second and third waveforms are less than 1 m. For the second waveform, the initial values are rough, while for the third waveform, RLD_A only provides the mixed peak position.
Figure 17 shows the fitting results. When the water depth is more than 1 m, the surface returns fitted by the two models are both correct, but the peak position of the surface return fitted by the EW model is closer to the reference data (see Figure 17a,d). When the water depth is less than 1 m, although the whole waveform is well fitted with the Gaussian model, the position of the surface return is not accurate (see Figure 17b,c). In contrast, the EW model is still applicable even in the absence of accurate initial values (see Figure 17e,f). In addition, it can be found that the fitting of water column scattering plays a very small role here, which also follows the practical situation, but it still influences the positions of the fitted surface return and bottom return.

3.3. Experiment Ⅲ: Performance Analysis for the Processing Algorithms in Deep Water

3.3.1. Reference Data

Dataset 2 is located on a slope with depths ranging from 20 m to 40 m, as shown in Figure 18. This dataset is selected to test the ability of the signal detection algorithms in deep water. The detection of surface signals has been analyzed in shallow water, and for the waveforms in deep water, the difficulty of signal detection mainly lies in the detection of the weak bottom signal. Thus, this experiment focuses on the detection of bottom signals.
Here, signals detected from the PMT3 channel are used as reference data. Because PMT3 has a larger field of view than PMT1, the intensity of the bottom return in PMT1 is weaker than that in PMT3 (see Figure 19b). In Dataset 2, the bottom signals are weak in PMT1, but still strong in PMT3. Therefore, we apply the proposed method to PMT1, and the bottom signals detected from PMT3 with the maximum detection method are used as reference data. However, waveforms in PMT3 are not applicable to all cases. When the water is shallow, PMT3 is affected by multiple reflection (see Figure 19a), and when the water depth reaches a certain value, the bottom return in PMT3 is also weak (see Figure 19c). Since the bottom reflection can be approximated as a non-directional diffuse scattering, there is no difference between the bottom return in PMT1 and PMT2. Therefore, the optimal strategy is to process the waveforms in PMT1 when the water is shallow and process the waveforms in PMT3 when the detection of the bottom signal in PMT1 has failed. However, it is still necessary to detect the weak bottom signal in PMT3 effectively. In Dataset 2, the bottom signals in PMT1 are weak, which can simulate the detection of the weak bottom signal in PMT3.

3.3.2. Bottom Signal Detection

In deep water, waveforms may have a weak bottom signal with an amplitude approximately equal to the noise level, which is possibly caused by a deep depth, a turbid water or a dark bottom. Therefore, besides the RLD, which is used to improve the signal resolution, a filtering algorithm which can de-noise the waveform but keep the bottom signal in the meantime is also effective. ASDF is used as a filtering algorithm in [23], which is a substitute for the direct cross-correlation function to estimate the time delay in two discrete time series and is more computationally convenient. ASDF_M and ASDF_A denote the point clouds generated from the signals detected by the maximum method with a fixed threshold and the adaptive threshold from the waveforms that are processed by ASDF. We also applied the method proposed in [25] which used a first derivative of a wide Gaussian filter to reduce the processing range of the waveform, and processed the waveforms in this range with normalization, multiple smoothing and three times derivation to increase the detection rate of the bottom signal. The point clouds generated from this method are referred to as dddNCFWF. For the waveform decomposition, since the Gaussian model is no longer applicable, the quadrilateral model (QUAD) introduced in [18] is applied here together with the EFSP model proposed in this paper. These models are initialized with RLD_A, and the corresponding point clouds are referred to as QUAD and EFSP. Before the performance analysis, all the point clouds are filtered according to neighboring points to eliminate outliers.
Figure 20 shows the detected water bottom on the selected profile (see Figure 18). The number of the detected points in ASDF_M is nearly twice that in Max or RLD_M. However, when the thresholds in the maximum method are adaptive, more water bottom points can be detected. RLD_A performs better than ASDF_A, which is exactly the opposite of RLD_M and ASDF_M. The intensity of the bottom signal varies greatly with water depth, but the fixed threshold cannot take all kinds of situations into account. The detection rate of ASDF_M is higher than Max and RLD_M, while the detection rate of ASDF_A is almost equal to that of RLD_A. Because noise can be filtered by ASDF, the results indicate that the fixed threshold is more sensitive to the noise than the adaptive threshold. dddNCFWF provides a high detection rate just as RLD_A. However, there is an offset between it and the reference data which may be induced by the dissimilarity between the emitted waveform and the symmetrical filters used in the method. Even with RLD_A as the initial value, QUAD still cannot accurately fit the waveforms with the depths over 35 m. In contrast, the EFSP model is more flexible, but it seems that there is little difference between RLD_A and EFSP.
In further comparing these algorithms, the above point clouds are statistically analyzed, as shown in Table 3; the detection rate (Dr) is the percentage of bottom points with an error less than sqrt(0.32+(0.015 × d)2) m (where d is the water depth); RMSE is the root mean square error of the bottom points; and max(d) denotes the detected maximum depth. Since there is an offset in dddNCFWF, its translation (dddNCFWF_T) that best matches the reference data is also considered in the following analysis.
The detection rates of Max and RLD_M are low, while ASDF_M shows a higher detection rate. The detection rate is significantly improved in ASDF_A, RLD_A and dddNCFWF. The detection rate of dddNCFWF is almost equal to ASDF_A but lower than RLD_A, because the fake signals with relatively high amplitude can be detected in dddNCFWF sometimes. In contrast, RLD can effectively keep the signal which is similar to the emitted waveform and remove other noise even with relatively high amplitude. The detection rate of QUAD is lower than RLD_A, indicating that the detection rate is decreased after waveform decomposition with the quadrilateral model. EFSP can improve the detection results in RLD_A with the highest detection rate of 76.64%. RMSE reflects the accuracy of the detected points. ASDF_M has the minimum RMSE, but its detection rate is too low. dddNCFWF has the largest RMSE. Although the accuracy can be improved by translation, there is still a large RMSE in dddNCFWF_T. The distance of the offset may be related to the signal strength, so it cannot be completely eliminated by the translation. Errors are introduced by the quadrilateral model with the RMSE increased from 0.1347 m to 0.1504 m in QUAD. After waveform decomposition with the EFSP model, the RMSE is reduced from 0.1347 m to 0.1076 m, indicating that EFSP may correct the initial value even if there is some error in it. EFSP has a great max(d), exceeding the traditional maximum method by 10 m. In addition, dddNCFWF has the maximum max(d) among the tested methods, implying that it may have a better performance with an appropriate filter.

3.3.3. Adaptability Analysis of the EFSP Model

Since water depth is a key factor affecting the strength of the bottom signal, four waveforms in different depths are selected to test the adaptability of the model in waveform decomposition. The amplitudes of the bottom returns in these waveforms are 21.4, 3.2, 2.5, and 2 bins, respectively. For the first three waveforms, the initial values provided by RLD_A are correct, while the position of the bottom signal in RLD_A is wrong for the last waveform.
The received waveforms and their fitting results are shown in Figure 21. It can be found that the quadrilateral model cannot fit water column scattering accurately compared with the EFSP model. When the bottom signal is fairly strong, the waveform decompositions with the two models can both converge to the exact position (see Figure 21a,b). With the decrease of bottom signal intensity, the inadaptability of the quadrilateral model becomes more and more obvious, even when the exact initial values are provided (see Figure 21c,e,g). In contrast, as long as the strengths of the bottom signals are sufficient, the locations of the bottom signals obtained with the EFSP model are precise regardless of the correctness of the initial values (see Figure 21d,f,h).

3.4. Experiment Ⅳ: Accuracy Assessment for the Processing Algorithms

3.4.1. Simulated Data

To evaluate the accuracy of the surface and bottom signals, the water LiDAR waveform model (Wa-LiD) presented by Abdallah et al. [42] was applied in this experiment. Wa-LiD is a successful simulator for simulating green channel waveforms received from water, which has been widely used in ALB research [15,18,19]. It can perfectly reproduce the received waveforms by adjusting some realistic water parameters [42]. In this experiment, the Gaussian function in Wa-LiD was replaced by the calibration waveform to better fit the real received waveforms. As shown in Figure 22, the real waveforms can be well fitted by Wa-LiD with proper environmental parameters. However, there are some differences between the leading and falling edges of the signal, which may be induced by the slope of water surface and bottom.
To assess the accuracy of the bathymetric signals both in shallow and deep water, we generated 10,000 simulated waveforms in which the depth varied from 0 to 2 m with an interval of 0.0002 m, and 10,000 waveforms in which the depth varied from 40 m to 50 m with an interval of 0.001 m. The environmental parameters of each simulated waveform were randomly selected within the normal range which was determined by the real waveforms.

3.4.2. Signal Detection in Shallow Water

For the simulated waveforms of shallow water (0–2 m), we test the waveform processing algorithms in Section 3.2.2. Since the accurate surface and bottom signal positions are known, the detection rate of surface (Dr_S) and bottom signal (Dr_B), and their RMSE (denoted by RMSE_S and RMSE_S, respectively) are assessed in Table 4. The accuracy of the surface signal is consistent with the field experiment results. EW still has the highest detection rate and accuracy. The Dr_S in simulated data is generally higher than that in field data, because the effect of waves on the water surface is not considered in the simulated waveforms. It is worth noting that the Dr_S of GD is higher than that of RLD_A, but the Dr_B of GD is significantly lower than that of RLD_A indicating that GD is not adaptive to the waveforms. Due to the changes of the environmental parameters, the Dr_S of RLD_A is lower than that of RLD_M. However, EW is still able to achieve the correct decomposition with inaccurate initial values, which depends on the appropriate model (EW) and a good optimization algorithm (TR).

3.4.3. Signal Detection in Deep Water

The waveform processing methods in Section 3.3.2 are evaluated with the simulated waveforms of deep water (40–50 m), and the statistical results are shown in Table 5. By comparing Table 3 and Table 5, it can be found that the translation of dddNCFWF (dddNCFWF_T) performs better than EFSP here. Because the simulated waveform does not take into account the signal stretching caused by the slope of the water surface, the offset in dddNCFWF becomes a fixed value which can be eliminated by translation. It also proves that the error of dddNCFWF is mostly caused by the mismatch between the symmetric filter and the asymmetric emitted waveform. Apart from dddNCFWF_T, EFSP still has the best performance. Consistent with the results of the field experiment, the adaptive threshold significantly increases the Dr_B, and the Dr_B of RLD_A is still greater than that of ASDF_A or dddNCFWF, showing that RLD_A are the best initial values. Although EFSP does not increase the Dr_B of RLD_A, the RMSE_B of EFSP has been significantly improved. Even for the weak bottom signals in deep water, constructing a suitable model is also useful for improving the accuracy.

4. Conclusions

In this paper, a depth-adaptive waveform decomposition method for ALB was developed by classifying and fitting the waveforms in the green channel based on the similarities between them and the water column scattering. The application of the proposed method in shallow water datasets (where the depth is less than 2 m) and deep water datasets (with depths ranging from 20 m to 60 m) showed that this method can cope with most of the waveforms and significantly improves the accuracy of the detected signals. The main conclusions are as follows:
  • Water column scattering can be used as a sign to distinguish the received waveforms in terms of depth. The defined parameter S can be used to measure the similarity between the received waveforms and the water column scattering. Since water column scattering is covered in the shallow water waveform, the S of the shallow water waveform is obviously greater than that of the deep water waveform. Thus, waveforms can be classified precisely according to S.
  • For the waveform preprocessing, improving the signal resolution is more efficient than denoising. With an appropriate signal detection threshold, RLD always performs better than ASDF with a higher signal detection rate. Although filtering algorithms can remove the noise in signals and improve the accuracy of signal detection, the weak bottom signal may be filtered out as noise in the meantime.
  • The adaptive threshold can improve the reliability of the signal detection. The intensity of the bottom signal varies greatly with water depth, while the noise in water column scattering may be stronger than the bottom signal, leading to the detection of fake signals. Furthermore, although RLD is a deconvolution algorithm with good noise resistance, noise is inevitably introduced in the process. The adaptive threshold can better cope with the fake signals because it takes into account the effects of the water column.
  • With an appropriate model and reliable initial values, waveform decomposition can significantly improve the signal detection rate and accuracy. The proposed models, EW and EFSP, can fit the waveforms well in most cases. Compared with the Gaussian function, the transformation of the calibration waveform can better fit the water surface and bottom signals. The exponential function with a second-order polynomial is consistent with the shape of water column scattering in the waveform. The TR algorithm can solve the model parameters in a reasonable region and provide an accurate solution. The results of waveform decomposition are based on the whole waveform and are accurate to the sub-sampling interval. Even when the initial values are wrong, the detection results can be corrected by waveform decomposition in some cases. In addition, the processing time of waveform decomposition is long, meaning that whether the wave decomposition step should be added depends on the accuracy requirements in practical applications. The waveform decomposition model proposed in this paper is for the Mapper5000 system and may need relevant adjustments when applied to waveforms acquired by other ALB systems.
The applicability of the model is very important in waveform decomposition. In shallow water, due to the random fitting of water column scattering, a bias between the surface return fitted by EW and the surface signal detected from the NIR channel can be found. For future research, it would be useful to study the fitting of water column scattering in extremely shallow water. Using physical models rather than empirical models to fit water column scattering may be a breakthrough in this field.

Author Contributions

Conceptualization, S.X., D.W., Q.X., Y.L. and P.L.; methodology, D.W., S.X., Q.X., Y.L. and P.L.; formal analysis, S.X., D.W., Q.X., Y.L., P.L., L.J., X.Z. and C.L.; investigation, S.X., D.W. and Q.X.; resources, S.X. and D.W.; data curation, S.X., L.J., X.Z. and C.L.; Writing—original draft preparation, S.X. and D.W.; Writing—review and editing, S.X., D.W., Q.X., Y.L., P.L., L.J., X.Z. and C.L.; visualization, D.W.; funding acquisition, S.X. and Q.X.

Funding

This research was funded by the National Natural Science Foundation of China (No. 41876105 and No. 41371436), the Foundation of Science and Technology on Near-surface Detection Laboratory (No. TCGZ2017A008), and the Information Engineering University Foundation for the Excellent Young (No. 2016610902).

Acknowledgments

The data used in this study were provided by Shanghai Institute of Fine Mechanics and Optics, Chinese Academy of Sciences. The positioning model of Mapper5000 was provided by Shandong University of Science and Technology. We are very grateful for their supports.

Conflicts of Interest

The authors declare no conflict of interest.

Appendix A

As shown in Figure A1, the NIR laser and the green laser share a scanner which produces an elliptical scan pattern, and the footprint has an almost elliptical shape at the sea surface [43,44]. With an altitude of 300 m, the footprint of the NIR laser is 0.8 m along-track and 0.78 m cross-track, and the footprint of the green laser is 0.32 m along-track and 0.31 m cross-track. The average point density on the sea bottom amounts to 0.8 points/m2. The IFOV of the receiver is divided into two parts: a small center IFOV (6 mrad) and a large edge IFOV (6–40 mrad). The APD, PMT1 and PMT2 channels use the small IFOV while the PMT3 channel uses the large IFOV. The sensor simultaneously collects full waveforms from four channels (NIR, PMT1–PMT3) with a sampling interval of 1 ns, corresponding to a distance of 0.15 m on the slant path in the air.
Figure A1. Schematic diagram of optical system in the Mapper5000 system [44]. IFOV: instantaneous field of view; APD: avalanche photodiode; PMT: photomultiplier tube; PBS: polarization beam splitter. The “Laser” contains a green laser and an NIR laser. The IFOV of the receiver is separated into two parts: a small center IFOV (6 mrad) and a large edge IFOV (6–40 mrad). The green light obtained from the small IFOV is split by a PBS and is received by PMT1 and PMT2.
Figure A1. Schematic diagram of optical system in the Mapper5000 system [44]. IFOV: instantaneous field of view; APD: avalanche photodiode; PMT: photomultiplier tube; PBS: polarization beam splitter. The “Laser” contains a green laser and an NIR laser. The IFOV of the receiver is separated into two parts: a small center IFOV (6 mrad) and a large edge IFOV (6–40 mrad). The green light obtained from the small IFOV is split by a PBS and is received by PMT1 and PMT2.
Sensors 19 05065 g0a1

Appendix B

As shown in Figure A2a, the scanner in the Mapper5000 system contains a rotating mirror which is controlled by a motor. The angle between the laser emission direction and the axis of rotation is 45°, and the angle between the normal of the reflecting mirror and the axis of rotation is 7.5°. As the mirror rotates, the laser will create an elliptical trajectory on the sea surface (see Figure A2b).
In this paper, the point clouds are generated by a data processing software which is designed for Mapper5000. After extracting the bathymetric signals from the waveform, a data file is generated containing the surface signal position tS, bottom signal position tB, the emitted signal position t0, the encoder data of the scanner and the GPS time of each waveform. The data file and the POS data are input into this software, and the point cloud data can be generated. The refraction correction is also completed in this software, the principle is as follows.
Figure A2. (a) The rotating scanner of Mapper5000 and (b) the spot trajectory on sea surface. The angle between the normal of the reflecting mirror and the axis of rotation is 7.5°, resulting in an elliptical spot trajectory on the sea surface [44]. (c) The mechanism of the refraction correction for Mapper5000. The ray path of the green laser in air and water are denoted by green and blue line, respectively.
Figure A2. (a) The rotating scanner of Mapper5000 and (b) the spot trajectory on sea surface. The angle between the normal of the reflecting mirror and the axis of rotation is 7.5°, resulting in an elliptical spot trajectory on the sea surface [44]. (c) The mechanism of the refraction correction for Mapper5000. The ray path of the green laser in air and water are denoted by green and blue line, respectively.
Sensors 19 05065 g0a2
As shown in Figure A2c, a right-handed coordinate system is established with the laser as the origin, the flight direction as the positive direction of the Y-axis, and vertical upward direction as the positive direction of the Z-axis. According to the encoder data, the emitted direction of the green laser can be determined, including the angle θa between the laser beam and the Z-axis, and the angle φ between the projection of the laser beam on the XY plane and the positive direction of the X-axis. Then the coordinates of the water surface point can be expressed as:
x S = O + r a r a = L a γ a = c a Δ t a 2 γ a = c a ( t S t 0 ) 2 γ a γ a = ( sin θ a cos ϕ , sin θ a sin ϕ , cos θ a ) T ,
where ra is the ray of the laser in the air equal to the product of its unit vector γa and the transmission distance La, ca is the velocity of light in air, and Δta is the two-way transmission time in the air. According to the Snell’s law, when the laser pulse incident on the air-water interface the incidence angle will change from θa to θw, and the velocity of light will decrease from ca to cw.
sin θ a sin θ w = n a n w = c w c a ,
where na and nw are the refractive indexes of air and water equal to 1 and 1.33, respectively. The coordinates of the water bottom point are given by:
x B = O + r a + r w r w = L w γ w = c w Δ t w 2 γ w = c w ( t B t S ) 2 γ w γ w = ( sin θ w cos ϕ , sin θ w sin ϕ , cos θ w ) T ,
where rw is the ray of the laser in the water equal to the product of its unit vector γw and the transmission distance Lw, and Δtw is the two-way transmission time in the water. Since a POS system is mounted on the same platform as the laser, the instantaneous position and attitude of the platform in WGS84 can be known from the GPS time and POS data, and the coordinates of the water surface and bottom points can be converted to WGS84.

References

  1. Hickman, G.D.; Hogg, J.E. Application of an airborne pulsed laser for near shore bathymetric measurements. Remote Sens. Environ. 1969, 1, 47–58. [Google Scholar] [CrossRef]
  2. Guenther, G.C. Airborne lidar bathymetry. In Digital Elevation Model Technologies and Applications: The Dem Users Manual, 2nd ed.; Maune, D.F., Ed.; ASPRS Publications: Bethesda, MD, USA, 2007; pp. 253–320. [Google Scholar]
  3. Guenther, G.C.; Cunningham, A.G.; Laroque, P.E.; Reid, D.J. Meeting the accuracy challenge in airborne Lidar bathymetry. In Proceedings of the 20th EARSeL Symposium: Workshop on Lidar Remote Sensing of Land and Sea, Dresden, Germany, 16–17 June 2000. [Google Scholar]
  4. Feurer, D.; Bailly, J.-S.; Puech, C.; Le Coarer, Y.; Viau, A.A. Very-high-resolution mapping of river-immersed topography by remote sensing. Prog. Phys. Geogr. 2008, 32, 403–419. [Google Scholar] [CrossRef]
  5. Collin, A.; Long, B.; Archambault, P. Salt-marsh characterization, zonation assessment and mapping through a dual-wavelength LiDAR. Remote Sens. Environ. 2010, 114, 520–530. [Google Scholar] [CrossRef]
  6. Finkl, C.W.; Makowski, C. Autoclassification versus cognitive interpretation of digital bathymetric data in terms of geomorphological features for seafloor characterization. J. Coast. Res. 2015, 31, 1–16. [Google Scholar] [CrossRef]
  7. Pan, Z.; Glennie, C.; Hartzell, P.; Fernandez-Diaz, J.; Legleiter, C.; Overstreet, B. Performance assessment of high resolution airborne full waveform lidar for shallow river bathymetry. Remote Sens. 2015, 7, 5133–5159. [Google Scholar] [CrossRef]
  8. Lesaignoux, A.; Bailly, J.; Feurer, D. Small water depth detection from green lidar simulated full waveforms: Application to gravel-bed river bathymetry. In Proceedings of the Physics in Signal and Image Processing (PSIP) Fifth Internation Conference, Mulhouse, France, 31 January–2 February 2007. [Google Scholar]
  9. Mandlburger, G.; Pfennigbauer, M.; Pfeifer, N. Analyzing near water surface penetration in laser bathymetry—A case study at the River Pielach. In Proceedings of the ISPRS Annals of the Photogrammetry, Remote Sensing and Spatial Information Sciences, Antalya, Turkey, 11–13 November 2013. [Google Scholar]
  10. Wang, C.-K.; Philpot, W.D. Using airborne bathymetric LiDAR to detect bottom type variation in shallow waters. Remote Sens. Environ. 2007, 106, 123–135. [Google Scholar] [CrossRef]
  11. Pe’eri, S.; Philpot, W. Increasing the existence of very shallow-water LIDAR measurements using the red-channel waveforms. IEEE Trans. Geosci. Remote Sens. 2007, 45, 1217–1223. [Google Scholar] [CrossRef]
  12. Allouis, T.; Bailly, J.S.; Pastol, Y.; Le Roux, C. Comparison of LiDAR waveform processing methods for very shallow water bathymetry using Raman, near-infrared and green signals. Earth Surf. Process. Landf. 2010, 35, 640–650. [Google Scholar] [CrossRef]
  13. Zhao, J.; Zhao, X.; Zhang, H.; Zhou, F. Shallow water measurements using a single green laser corrected by building a near water surface penetration model. Remote Sens. 2017, 9, 426. [Google Scholar] [CrossRef]
  14. Mallet, C.; Bretar, F. Full-waveform topographic lidar: State-of-the-art. ISPRS J. Photogramm. Remote Sens. 2009, 64, 1–16. [Google Scholar] [CrossRef]
  15. Wang, C.; Li, Q.; Liu, Y.; Wu, G.; Liu, P.; Ding, X. A comparison of waveform processing algorithms for single-wavelength LiDAR bathymetry. ISPRS J. Photogramm. Remote Sens. 2015, 101, 22–35. [Google Scholar] [CrossRef]
  16. Collin, A.; Archambault, P.; Long, B. Mapping the shallow water seabed habitat with the SHOALS. IEEE Trans. Geosci. Remote Sens. 2008, 46, 2947–2955. [Google Scholar] [CrossRef]
  17. Abdallah, H.; Bailly, J.S.; Baghdadi, N.; Saint-Geours, N.; Fabre, F. Potential of space-borne LiDAR sensors for global bathymetry in coastal and inland waters. IEEE J. Sel. Top. Appl. Earth Obs. Remote Sens. 2013, 6, 202–216. [Google Scholar] [CrossRef]
  18. Abady, L.; Bailly, J.S.; Baghdadi, N.; Pastol, Y.; Abdallah, H. Assessment of quadrilateral fitting of the water column contribution in lidar waveforms on bathymetry estimates. IEEE Geosci. Remote Sens. 2014, 11, 813–817. [Google Scholar] [CrossRef]
  19. Ding, K.; Li, Q.; Zhu, J.; Wang, C.; Guan, M.; Chen, Z.; Yang, C.; Cui, Y.; Liao, J. An improved quadrilateral fitting algorithm for the water column contribution in airborne bathymetric lidar waveforms. Sensors 2018, 18, 552. [Google Scholar] [CrossRef]
  20. Schwarz, R.; Mandlburger, G.; Pfennigbauer, M.; Pfeifer, N. Design and evaluation of a full-wave surface and bottom-detection algorithm for LiDAR bathymetry of very shallow waters. ISPRS J. Photogramm. Remote Sens. 2019, 150, 1–10. [Google Scholar] [CrossRef]
  21. Saylam, K.; Brown, R.A.; Hupp, J.R. Assessment of depth and turbidity with airborne Lidar bathymetry and multiband satellite imagery in shallow water bodies of the Alaskan North Slope. Int. J. Appl. Earth Obs. Geoinf. 2017, 58, 191–200. [Google Scholar] [CrossRef]
  22. Wu, J.; Van Aardt, J.; Asner, G.P. A comparison of signal deconvolution algorithms based on small-footprint lidar waveform simulation. IEEE Trans. Geosci. Remote Sens. 2011, 49, 2402–2414. [Google Scholar] [CrossRef]
  23. Wagner, W.; Roncat, A.; Melzer, T.; Ullrich, A. Waveform analysis techniques in airborne laser scanning. In Proceedings of the ISPRS Workshop on Laser Scanning 2007 and SilviLaser 2007, Espoo, Finland, 12–14 September 2007. [Google Scholar]
  24. Richter, K.; Maas, H.-G.; Westfeld, P.; Weiß, R. An Approach to Determining Turbidity and Correcting for Signal Attenuation in Airborne Lidar Bathymetry. PFG J. Photogramm. Remote Sens. Geoinf. Sci. 2017, 85, 31–40. [Google Scholar] [CrossRef]
  25. Launeau, P.; Giraud, M.; Robin, M.; Baltzer, A. Full-Waveform LiDAR Fast Analysis of a Moderately Turbid Bay in Western France. Remote Sens. 2019, 11, 117. [Google Scholar] [CrossRef]
  26. Parrish, C.E.; Jeong, I.; Nowak, R.D.; Brent Smith, R. Empirical comparison of full-waveform lidar algorithms: Range extraction and discrimination performance. Photogramm. Eng. Remote Sens. 2011, 77, 825–838. [Google Scholar] [CrossRef]
  27. Jutzi, B.; Stilla, U. Range determination with waveform recording laser systems using a Wiener Filter. ISPRS J. Photogramm. Remote Sens. 2006, 61, 95–107. [Google Scholar] [CrossRef]
  28. Richardson, W.H. Bayesian-based iterative method of image restoration. J. Opt. Soc. Am. 1972, 62, 55–59. [Google Scholar] [CrossRef]
  29. Lucy, L.B. An iterative technique for the rectification of observed distribution. Astron. J. 1974, 79, 745–775. [Google Scholar] [CrossRef]
  30. Wu, J.; Van Aardt, J.; McGlinchy, J.; Asner, G.P. A robust signal preprocessing chain for small-footprint waveform lidar. IEEE Trans. Geosci. Remote Sens. 2012, 50, 3242–3255. [Google Scholar] [CrossRef]
  31. Naik, R.K.; Sahu, P.K. Comprehensive study on deconvolution and denoising of LiDAR back-scattered signal. In Proceedings of the 2013 International Conference on Microwave and Photonics (ICMAP), Dhanbad, India, 13–15 December 2013. [Google Scholar]
  32. White, R.L. Image restoration using the damped Richardson-Lucy method. In The Restoration of HST Images and Spectra-II; Hanisch, R.J., White, R.L., Eds.; STScI: Baltimore, MD, USA, 1994. [Google Scholar]
  33. Wagner, W.; Ullrich, A.; Melzer, T.; Briese, C.; Kraus, K. From single-pulse to full-waveform airborne laser scanners: Potential and practical challenges. In Proceedings of the International Archives of Photogrammetry, Remote Sensing and Spatial Information Sciences, Istanbul, Turkey, 12–23 July 2004; pp. 201–206. [Google Scholar]
  34. Schwarz, R.; Pfeifer, N.; Pfennigbauer, M.; Ullrich, A. Exponential decomposition with implicit deconvolution of lidar backscatter from the water column. PFG J. Photogramm. Remote Sens. Geoinf. Sci. 2017, 85, 159–167. [Google Scholar] [CrossRef]
  35. Persson, Å.; Söderman, U.; Töpel, J.; Ahlberg, S. Visualization and analysis of full-waveform airborne laser scanner data. Int. Arch. Photogramm. Remote Sens. Spat. Inf. Sci. 2005, 36, 103–108. [Google Scholar]
  36. Hofton, M.A.; Minster, J.B.; Blair, J.B. Decomposition of laser altimeter waveforms. IEEE Trans. Geosci. Remote Sens. 2000, 38, 1989–1996. [Google Scholar] [CrossRef]
  37. Li, P.C.; Xu, Q.; Xing, S.; Liu, Z.; Geng, X.; Hou, X.; Zhang, J. Full-waveform LiDAR data decomposition method based on global convergent LM. Infrared Laser Eng. 2015, 44, 2262–2267. [Google Scholar]
  38. Hernandez-Marin, S.; Wallace, A.M.; Gibson, G.J. Bayesian analysis of lidar signals with multiple returns. IEEE Trans. Pattern Anal. Mach. Intell. 2007, 29, 2170–2180. [Google Scholar] [CrossRef]
  39. Yuan, Y. Recent advances in trust region algorithms. Math. Program. 2015, 151, 249–281. [Google Scholar] [CrossRef]
  40. Wang, X.; Yuan, Y. An augmented Lagrangian trust region method for equality constrained optimization. Optim. Method Softw. 2015, 30, 559–582. [Google Scholar] [CrossRef]
  41. Coleman, T.F.; Li, Y. An interior trust region approach for nonlinear minimization subject to bounds. SIAM J. Control. 1996, 6, 418–445. [Google Scholar] [CrossRef] [Green Version]
  42. Abdallah, H.; Baghdadi, N.; Bailly, J.S.; Pastol, Y.; Fabre, F. Wa-LiD: A new LiDAR simulator for waters. IEEE Geosci. Remote Sens. Lett. 2012, 9, 744–748. [Google Scholar] [CrossRef] [Green Version]
  43. Hu, S.; He, Y.; Chen, W.; Zhu, X.; Zang, H.; Lyu, D.; Tian, M.; Yu, J.; Tao, B.; Huang, T.; et al. Design of airborne dual-frequency laser radar system. Infrared Laser Eng. 2018, 47, 0930001. [Google Scholar]
  44. He, Y.; Hu, S.; Chen, W.; Zhu, X.; Wang, Y.; Yang, Z.; Zhu, X.; Lyu, D.; Yu, J.; Huang, T.; et al. Research Progress of Domestic Airborne Dual-Frequency LiDAR Detection Technology. Laser Optoelectron. Prog. 2018, 55, 082801. [Google Scholar]
Figure 1. Study area. The yellow solid lines denote the distribution of the strips, the range of the datasets are in the red rectangles, and “+” denote the positions of the bathymetric points (P1–P4).
Figure 1. Study area. The yellow solid lines denote the distribution of the strips, the range of the datasets are in the red rectangles, and “+” denote the positions of the bathymetric points (P1–P4).
Sensors 19 05065 g001
Figure 2. Workflow of the depth-adaptive waveform decomposition method. PMT: photomultiplier tube; RLD: Richardson–Lucy deconvolution; EFSP: exponential function with second-order polynomial model; EW: empirical waveform model; NIR: near-infrared; TR: trust region algorithm.
Figure 2. Workflow of the depth-adaptive waveform decomposition method. PMT: photomultiplier tube; RLD: Richardson–Lucy deconvolution; EFSP: exponential function with second-order polynomial model; EW: empirical waveform model; NIR: near-infrared; TR: trust region algorithm.
Sensors 19 05065 g002
Figure 3. The useful range of the received waveform. The received waveform recorded by the ALB system ranges from 0–3200 ns, while the useful range is 2389–2562 ns.
Figure 3. The useful range of the received waveform. The received waveform recorded by the ALB system ranges from 0–3200 ns, while the useful range is 2389–2562 ns.
Sensors 19 05065 g003
Figure 4. The received waveforms with depths of (a) 0.6 m, (b) 10.8 m and (c) 23.9 m.
Figure 4. The received waveforms with depths of (a) 0.6 m, (b) 10.8 m and (c) 23.9 m.
Sensors 19 05065 g004
Figure 5. R(t) of (a) a shallow water (SW) waveform with a large S of 4220.5 and (b) a deep water (DW) waveform with a small S of 133.4.
Figure 5. R(t) of (a) a shallow water (SW) waveform with a large S of 4220.5 and (b) a deep water (DW) waveform with a small S of 133.4.
Sensors 19 05065 g005
Figure 6. Examples of the deconvoluted waveforms by RLD in (a) shallow water and (b) deep water. RLD: Richardson–Lucy Deconvolution. The result of the deep water waveform is only displayed in the portion near the bottom signal.
Figure 6. Examples of the deconvoluted waveforms by RLD in (a) shallow water and (b) deep water. RLD: Richardson–Lucy Deconvolution. The result of the deep water waveform is only displayed in the portion near the bottom signal.
Sensors 19 05065 g006
Figure 7. Comparison of the conventional signal detection methods, where “threshold1” used a threshold that was equal to one-tenth of the maximum amplitude of wR(t) (MA), “threshold2” used a threshold that was equal to 0.05 MA, “maximum1” used a threshold that was equal to 0.2 MA, and “maximum2” used a threshold that was equal to 0.1 MA.
Figure 7. Comparison of the conventional signal detection methods, where “threshold1” used a threshold that was equal to one-tenth of the maximum amplitude of wR(t) (MA), “threshold2” used a threshold that was equal to 0.05 MA, “maximum1” used a threshold that was equal to 0.2 MA, and “maximum2” used a threshold that was equal to 0.1 MA.
Sensors 19 05065 g007
Figure 8. Thresholds used in signal detection. The high-intensity fake signals are located in water column scattering, which may be detected by the maximum method with a fixed threshold.
Figure 8. Thresholds used in signal detection. The high-intensity fake signals are located in water column scattering, which may be detected by the maximum method with a fixed threshold.
Sensors 19 05065 g008
Figure 9. Calibration waveform (solid) and fitting curves (dashed) with Gaussian, exponential and harmonic models. The Gaussian model includes only one Gaussian function. The exponential model is a sum of two exponential functions, and the harmonic model is a combination of four harmonic functions. All fitting parameters were solved by a Levenberg–Marquardt algorithm, and the adjusted R-Square (R2) coefficients of Gaussian, exponential and harmonic models are 0.9332, 0.8015 and 0.9956 respectively.
Figure 9. Calibration waveform (solid) and fitting curves (dashed) with Gaussian, exponential and harmonic models. The Gaussian model includes only one Gaussian function. The exponential model is a sum of two exponential functions, and the harmonic model is a combination of four harmonic functions. All fitting parameters were solved by a Levenberg–Marquardt algorithm, and the adjusted R-Square (R2) coefficients of Gaussian, exponential and harmonic models are 0.9332, 0.8015 and 0.9956 respectively.
Sensors 19 05065 g009
Figure 10. Water column scattering model.
Figure 10. Water column scattering model.
Sensors 19 05065 g010
Figure 11. (a) Mean curves and (b) variance curves of the four bathymetric points. The time offset indicates the offset of the waveform from the peak of the surface signal.
Figure 11. (a) Mean curves and (b) variance curves of the four bathymetric points. The time offset indicates the offset of the waveform from the peak of the surface signal.
Sensors 19 05065 g011
Figure 12. The distribution of S in different cases.
Figure 12. The distribution of S in different cases.
Sensors 19 05065 g012
Figure 13. The scope of Dataset 1 (red line) and the selected profile (yellow line).
Figure 13. The scope of Dataset 1 (red line) and the selected profile (yellow line).
Sensors 19 05065 g013
Figure 14. The (a) normal, (b) buried and (c) saturated NIR signals. The NIR signals are not always reliable. Sometimes the surface signal is buried below the noise level which cannot be detected, because its strength is not within the dynamic range of the receiver. Sometimes the surface signal is saturated, because the receiving direction of the sensor is approximately in the reflection direction of sunlight.
Figure 14. The (a) normal, (b) buried and (c) saturated NIR signals. The NIR signals are not always reliable. Sometimes the surface signal is buried below the noise level which cannot be detected, because its strength is not within the dynamic range of the receiver. Sometimes the surface signal is saturated, because the receiving direction of the sensor is approximately in the reflection direction of sunlight.
Sensors 19 05065 g014
Figure 15. Comparison between the height of the water surface in reference data and detection results.
Figure 15. Comparison between the height of the water surface in reference data and detection results.
Sensors 19 05065 g015
Figure 16. Locations of the water surface on the selected profile from: (a) reference data; (b) Max; (c) RLD_M; (d) RLD_A; (e) Gaussian decomposition (GD); and (f) EW.
Figure 16. Locations of the water surface on the selected profile from: (a) reference data; (b) Max; (c) RLD_M; (d) RLD_A; (e) Gaussian decomposition (GD); and (f) EW.
Sensors 19 05065 g016
Figure 17. Waveforms fitted with (ac) the Gaussian model and (df) the EW model.
Figure 17. Waveforms fitted with (ac) the Gaussian model and (df) the EW model.
Sensors 19 05065 g017
Figure 18. The scope of Dataset 2 (red line) and the selected profile (yellow line).
Figure 18. The scope of Dataset 2 (red line) and the selected profile (yellow line).
Sensors 19 05065 g018
Figure 19. The received waveforms in PMT1, PMT2 and PMT3 with depths of (a) 1.4 m, (b) 35.1 m and (c) 57.8 m.
Figure 19. The received waveforms in PMT1, PMT2 and PMT3 with depths of (a) 1.4 m, (b) 35.1 m and (c) 57.8 m.
Sensors 19 05065 g019
Figure 20. Locations of the water bottom on the selected profile from: (a) Max; (b) ASDF_M; (c) RLD_M; (d) ASDF_A (e) RLD_A; (f) dddNCFWF; (g) quadrilateral model (QUAD) and (h) EFSP.
Figure 20. Locations of the water bottom on the selected profile from: (a) Max; (b) ASDF_M; (c) RLD_M; (d) ASDF_A (e) RLD_A; (f) dddNCFWF; (g) quadrilateral model (QUAD) and (h) EFSP.
Sensors 19 05065 g020
Figure 21. Waveforms fitted with (a,c,e,g) the quadrilateral model and (b,d,f,h) the EFSP model (zoom in the bottom). The peak of the fitted bottom return is marked by “+”.
Figure 21. Waveforms fitted with (a,c,e,g) the quadrilateral model and (b,d,f,h) the EFSP model (zoom in the bottom). The peak of the fitted bottom return is marked by “+”.
Sensors 19 05065 g021
Figure 22. Comparison between the real waveforms and the simulated waveforms with depths of (a) 0.8 m and (b) 27 m.
Figure 22. Comparison between the real waveforms and the simulated waveforms with depths of (a) 0.8 m and (b) 27 m.
Sensors 19 05065 g022
Table 1. Field data acquisition parameters.
Table 1. Field data acquisition parameters.
Pulse Repetition FrequencySampling SpeedBeam DivergenceAltitudeSpeedSwath Width
5 kHz1 GHzNIR: 2.5 mrad
green: 1 mrad
300 m190 km/h160 m
Table 2. Performance assessments for surface detection results with field data.
Table 2. Performance assessments for surface detection results with field data.
AlgorithmDr (%)RMSE (m)min(d) (m)Std.
Max95.390.16450.44140.2276
RLD_M97.530.12180.32920.1309
RLD_A97.190.10850.32920.1107
GD97.290.11740.22080.1456
EW99.110.09010.22080.1159
Reference58.540.0839
Table 3. Performance assessments for bottom detection results with field data.
Table 3. Performance assessments for bottom detection results with field data.
AlgorithmDr (%)RMSE (m)max(d) (m)
Max31.220.169429.57
ASDF_M53.800.094435.25
RLD_M33.130.149830.36
ASDF_A71.160.112437.50
RLD_A73.880.134738.60
dddNCFWF71.010.379841.14
dddNCFWF_T71.980.224740.95
QUAD66.390.150437.36
EFSP74.640.107640.49
Table 4. Performance assessments for detection results in shallow water with simulated data.
Table 4. Performance assessments for detection results in shallow water with simulated data.
AlgorithmDr_S (%)Dr_B (%)RMSE_S (m)RMSE_B (m)min(d) (m)
Max69.8866.960.14330.63650.5920
RLD_M76.8174.950.14370.36860.4178
RLD_A71.7774.870.13110.12660.4180
GD78.7059.070.16970.88320.1540
EW94.7597.920.10590.08450.0558
Table 5. Performance assessments for detection results in deep water with simulated data.
Table 5. Performance assessments for detection results in deep water with simulated data.
AlgorithmDr_S (%)Dr_B (%)RMSE_S (m)RMSE_B (m)max(d) (m)
Max1000.000.0469
ASDF_M1000.220.05080.083543.36
RLD_M10020.350.08570.189649.92
ASDF_A10042.290.05080.101249.70
RLD_A10056.340.08570.170249.92
dddNCFWF10052.290.39940.421749.92
dddNCFWF_T10059.610.04340.085049.92
QUAD10046.970.14160.257349.92
EFSP10056.690.06160.068149.92

Share and Cite

MDPI and ACS Style

Xing, S.; Wang, D.; Xu, Q.; Lin, Y.; Li, P.; Jiao, L.; Zhang, X.; Liu, C. A Depth-Adaptive Waveform Decomposition Method for Airborne LiDAR Bathymetry. Sensors 2019, 19, 5065. https://doi.org/10.3390/s19235065

AMA Style

Xing S, Wang D, Xu Q, Lin Y, Li P, Jiao L, Zhang X, Liu C. A Depth-Adaptive Waveform Decomposition Method for Airborne LiDAR Bathymetry. Sensors. 2019; 19(23):5065. https://doi.org/10.3390/s19235065

Chicago/Turabian Style

Xing, Shuai, Dandi Wang, Qing Xu, Yuzhun Lin, Pengcheng Li, Lin Jiao, Xinlei Zhang, and Chenbo Liu. 2019. "A Depth-Adaptive Waveform Decomposition Method for Airborne LiDAR Bathymetry" Sensors 19, no. 23: 5065. https://doi.org/10.3390/s19235065

APA Style

Xing, S., Wang, D., Xu, Q., Lin, Y., Li, P., Jiao, L., Zhang, X., & Liu, C. (2019). A Depth-Adaptive Waveform Decomposition Method for Airborne LiDAR Bathymetry. Sensors, 19(23), 5065. https://doi.org/10.3390/s19235065

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop