Open Access
Add to BibSonomy
Abstract
Determination of the correct water surface height (WSH) from green laser (532 nm) echoes alone in bathymetry LiDAR is challenging, as the green laser return near the water surface involves both specular reflection from the air–water interface and backscattered return from the water volume. In this paper, a low-complexity method based on linear approximation of the leading edge (LLE) is proposed. The results of this LLE method were compared with those of three common algorithms of peak detection, half peak power, and surface-volume-bottom implemented on airborne datasets with various surface roughness conditions. In addition, the method was evaluated in waters with a wide range of optical properties through a controllable tank experiment. The uncertainties in the WSHs of all algorithms were greater when the water volume backscattering dominated the surface return; they were sensitive to variations in the optical properties of water, and increased exponentially with decreasing LiDAR attenuation coefficient (KLiDAR). Comparatively, the LLE algorithm had the fastest computational speed and demonstrated the best performance in situations where specular reflection or volume backscatter return was dominant, with average and maximum errors of less than 0.06 and 0.13 m, respectively.
© 2022 Optica Publishing Group under the terms of the Optica Open Access Publishing Agreement
1. Introduction
Airborne LiDAR bathymetry (ALB) is a proven technology for measuring the depths of relatively shallow coastal waters from the air by using a scanning, pulsed laser beam [1]. ALB uses a short green pulsed laser (typically 532 nm), which is capable of penetrating the water to measure the depth of a shallow coastal area from a low-altitude aircraft. The depth measurement is generally based on the round-trip transit time interval between the two echoes from the water surface and the bottom. The exact water surface height (WSH) or range in air between the LiDAR and water surface is required as the basis for measuring the water depth correctly, which is an inevitable precondition for performing the proper range and refraction correction of the laser pulse passing through two media [2].
However, owing to the capability of the green signal to penetrate water, the reflection of the green laser signal near the water surface involves specular reflection from the air–water interface and a diffused, backscattered return from the water column immediately below [2][3], which is commonly termed as “surface return”. The peak position of the LiDAR-detected surface return is an important feature for surface detection, but it usually cannot represent the location of the water surface exactly, resulting in an apparent water surface location (AWSL). The interaction of the laser beam with the air-water interface is complex, wherein the AWSL is influenced by various parameters, particularly the inherent optical properties (IOPs) of the water body, incidence angle, and water surface roughness [2]. Additionally, if the water surface is perfectly flat and mirror-like, which is more favorable for LiDAR bathymetry survey, weak or even no interface return would reach the transceiver because of the off-nadir geometry. Thus, Guenther theoretically defined a ratio between the specular reflection power (Ps) from the air–water interface and volume backscattering power (Pc) from the water column to determine the characteristics of the water surface return [3]. When the ratio of Ps/Pc is very small, that is, interface return is weaker than the water volume return, the WSH is often lower than the actual elevation when derived from the green signal alone based on its peak position, which directly causes significant errors in the derived water depths and laser spot positions at the water bottom. Guenther’s analysis suggested that the WSH error for reasonably clean water is typically in excess of 30 cm and can be as large as 60 cm; moreover, this error will be greater in the case of cleaner water and can result in significant uncertainties in depth retrieval [3]. The LiDAR system in a previous study [4] used the primary near-infrared (NIR; 1064 nm) pulse to detect the surface, as the depth error between the interface and volume returns for this design is smaller owing to reduced penetration of the water column at this wavelength. However, current small-foot-print bathymetry LiDARs no longer use the 1064 nm signal but emit and receive only the 532 nm green signal. Further, the slanting beam angle of ca. 20° usually employed in LiDAR bathymetry will not create adequate backscattering of the energy to the receiver to enable detection of the water surface echo for each 1064 nm pulse owing to reduced surface return probability [5], particularly during low wind conditions. Hence, it is necessary to determine the correct WSHs from the green laser echoes alone.
Mandlburger et al. conducted experiments at the Pielach River in Austria by using RIEGL VQ-820-G (532 nm) and VQ-580 (1064 nm) scanners. They found that the mean penetration of the green laser signal into the water column was approximately 0.25 m with the NIR signal as reference, and confirmed that the WSH is often underestimated when derived from the laser echoes alone [2]. It was also suggested that in the derivation of the water surface, the upper hull of the green laser echoes coincides with the NIR signal and restricts the underestimation of the water level to below 6 cm. However, these results should be obtained through statistical analysis of the aggregated neighboring echoes, which is feasible only for stationary water bodies such as ponds near the Pielach River. The actual position of the water surface may be very close to the leading edge of the observed waveforms based on the investigation in [6]. Pan et al. [7] proposed a simple leading-edge detection method based on a single threshold to accurately estimate the actual water surface, but its underestimation of the WSH can be up to 0.6 m in some regions, which was considered to be associated with water turbidity. Furthermore, Schwarz et al. [6] indicated that the lower the threshold, the better will be the concurrence between the estimated and actual positions of the water surface, but the noise rejection will be poor, as only a small fraction of the signal energy is utilized [8]. Later, Schwarz et al. [9] presented a surface-volume-bottom (SVB) algorithm, which accurately regarded the point of onset of the backscatter cross section of the water column as the position of the water surface. The SVB method can achieve smoother results when compared with the single threshold method, but it is more time-consuming, and the non-uniqueness of the numerical solution cannot be completely avoided. Thus, the feasibility of deriving the water surface from the green signal alone remains challenging; more importantly, the influence of various rough surface conditions and optical properties of water on the WSH determination requires extensive quantitative evaluation.
Therefore, in the present study, we conducted a series of airborne LiDAR tests and a controllable tank experiment to examine the performances of various algorithms (peak detection (PD), half peak power (HPP), and SVB)) and reveal their WSH uncertainties when the derivation is performed using green laser echoes in water bodies with various rough surface conditions and optical properties of water. Furthermore, we proposed a low-complexity but efficient method based on linear approximation of the leading edge (LLE) to accurately determine the WSH from green laser echoes alone.
The remainder of this study is organized as follows. Section 2 describes the airborne LiDAR data used in this study and the experimental approach. Section 3 presents a comparative evaluation of various algorithms used for water surface detection. Section 4 discusses the experimental results. Finally, Section 6 presents the concluding remarks.
2. Data and methods
2.1 Airborne LiDAR data
Mapper5000 shown in Fig. 1 is an airborne LiDAR system designed and manufactured by Shanghai Institute of Optics and Fine Mechanics (SIOM) [10][11]. The LiDAR system is mainly composed of three parts: oceanic LiDAR, cooling device, and integrated navigation unit (NovAtel SPAN-ISA-100C). Mapper5000 is a dual-wavelength LiDAR system, which simultaneously emits a 532 nm green laser beam and a coaxial 1064 nm NIR laser beam. The energies of the green and NIR laser beams are 3 mJ and 1 mJ, respectively, and the pulse repetition rate is 5 kHz. The scanner of Mapper5000 has a nutating design with the mirror axis located at an approximate offset of 7.5° from the axis of rotation (Fig. 1(b)). This scan mechanism produces a variable incident angle of the output laser beam with respect to the nadir from 10.5° to 15°, and creates a pseudo-elliptical scanning pattern on the ground. Mapper5000 employs segmented field-of-view (FOV) receivers, the optical design of which was first introduced in the CZMIL system [12][13] to detect the return signals over a wide dynamic range. Briefly, the return signals are focused by a lens onto a field separator, which divides the full FOV (40 mrad) into two parts comprising a small FOV (6 mrad) at the center and a ring-shaped FOV (6–40 mrad). Although the incident laser beam is scattered in the water column and spreads out spatially as an expanding cone whose size increase exponentially with increase in the optical depth [14][15], its energy on the water surface and at a small physical depth or optical thickness is mostly distributed across the small FOV channel. Thus, the waveforms from the small FOV channel collect most of the green laser return signals from the water surface and subsurface water column, which are used primarily for the analysis of the WSH. In addition, there is another channel for detecting the NIR returns. It is universally known that the penetration of the NIR laser pulse into the water column is limited to a few millimeters owing to the high absorption of water at 1064 nm; hence, only the direct specular reflections from the air-water interface can be detected. Because the receiving optical paths of the NIR and green channels coincide, the water surface returns from the NIR channel can then be used to derive the actual water surface level, which can be regarded as a reference for the evaluation of the results derived from the green small FOV channel.
Fig. 1. (a) Mapper5000 system comprising the cooling device, oceanic LiDAR, and integrated navigation system. (b) Schematic of nutating scanner in Mapper5000.
Airborne LiDAR measurements using Mapper5000 were conducted in three coastal areas of the East Island (109°24′50″E, 18°13′12″N), Dazhou Island (110°23′28″E, 18°47′43″N), and Ganquan Island (111°35′10″E, 16°30′28″N) of Hainan Province in August 2020 and October 2021. Table 1 lists the characteristics of the airborne Mapper5000 LiDAR bathymetry datasets. Partial airborne LiDAR data in the areas (shown as red rectangles in Fig. 2) several kilometers away from the three islands were selected to study the interaction of the laser pulses with the water surface. In these areas, the water generally has sufficient depth, and the water surface during both experiments was considered relatively calm but with some roughness, which made it possible to determine the water surface level more accurately with the NIR returns as reference. The water in the selected experimental areas was considered optically homogeneous. The LiDAR attenuation coefficients (KLiDAR) of the experimental areas near East Island, Dazhou Island, and Ganquan Island were approximately 0.29 m−1, 0.19 m−1, and 0.08 m−1, respectively.
Fig. 2. Experimental sites of (a) East Island, (b) Dazhou Island, and (c) Ganquan Island in Hainan Province of China. The flight trajectories are shown in yellow lines, and the selected experimental areas are denoted as red rectangles.
Table 1. Characteristics of airborne Mapper5000 LiDAR bathymetry experiment datasets
2.2 Tank experiments
To assess the influence of the optical property of water on the WSH derivation, control experiments were conducted in a tank with water of varying optical properties. A shipborne oceanic LiDAR (Fig. 3(a)) manufactured by SIOM was used in these experiments [16]. This LiDAR used the same optical filter, photomultiplier tube (PMT) detectors, and digitizer electronic subsystem as Mapper5000, but the major differences between the two devices were the energy, repetition rate of the 532 nm laser pulse, and receiving apertures (see Table 2). Briefly, a 532 nm linearly polarized laser with repetition rate of 100 Hz and energy per pulse of 50 µJ was used in the transmitting system, and an adjustable aperture was driven by a stepper motor to set the receiving FOV from 32 mrad to 105 mrad. Considering that the system impulse response functions of the two LiDAR were similar because their laser devices were manufactured by the same manufacturer, the waveforms from the shipborne oceanic LiDAR can fully represent the characteristics of those from Mapper5000. Moreover, a custom water tank with the size of 1.2 W × 1.2 D × 1.5 H m was filled with water to a depth of 1.5 m. To minimize the bottom of the 532 nm laser pulse reflection and stray light reflected from the tank walls, the inner and outer surfaces of the water tank were dull polished and covered with low-reflectivity black coating. An in-situ absorption and attenuation meter (AC-S) manufactured by WETLabs was placed in this tank (Fig. 3(b)) to measure the absorption coefficient (a) and beam attenuation coefficient (c) of the water body over a spectral range of 400–740 nm with 3.5 nm spectral interval. The value of c measured by the AC-S was corrected by adding the attenuation coefficient of pure water (0.0467 m−1 at 532 nm).
Fig. 3. (a) Device diagram of LiDAR maintaining the 15° off-nadir angle and loaded on the metal frame. (b) Tank experiment with an AC-S at the corner of the tank. (c-d) Schematics of different measurement methods.
Table 2. System parameters of Mapper5000 and shipborne LiDAR
The schematics in Fig. 3(c)-(d) illustrate the experimental processes. During the tank experiments, the shipborne LiDAR was fixed on a 16-m-high platform and the green laser pulse emitted from the LiDAR hit the water surface of the tank at an off-nadir incident angle of 15°, which was set to match the setting of the airborne LiDAR. Because the water surface during the tank experiments was mostly calm as a mirror, the specular reflection from the air–water interface could be avoided for ensuring green “surface” return dominated by the backscattered signals from the water column. The magnitude of interface specular reflection return was generally reduced below the detectable level. Additionally, there is no NIR channel in the shipborne LiDAR to derive the WSH for the reference. To obtain the exact water surface location as the reference WSH, a plastic plate (only 5 mm thick) was placed on the water surface. Owing to its small thickness, the surface of the plate was considered to be at the same height as the water surface. Note that the plastic plate was opaque to avoid the influence of the water column, and its surface was covered by dark gray paper to restrict the reflection intensity of the plate to within the dynamic range of the PMT detector and avoid signal saturation, the peak position of the measured plate waveform can be considered as the reference for the water surface, which can be found in Section 4.2.
More importantly, the IOPs (a and c) of the water in the tank were altered by adding different amounts of a sediment. After each addition of the sediment, the water was mixed well by the water pump on the bottom of the tank. Then, we turned off the water pump and conducted the LiDAR and AC-S measurements after the water surface became calm. During each measurement, we first recorded the waveforms from the water column to avoid the influence of sediment resettlement on the measurement to the extent possible; then, we moved the floating plastic plate to the laser spot for the acquisition of the reference WSH. Finally, these waveforms from the plastic plate and water column as well as the synchronous IOP measurements were used to quantitatively analyze the differences in the relationship between the WSL and IOPs. The detailed results of these experiments are presented and discussed in Section 4.2.
2.3 Estimation of Ps/Pc ratio
The ratio between Ps and Pc was theoretically defined by Guenther [3]. In this study, we proposed a practical method to roughly estimate Ps/Pc from the Mapper5000 waveforms. This method is based on the following two assumptions: (1) the laser attenuation process in the water body can be approximately expressed as an exponential function; (2) Ps/Pc can be approximated as the ratio of values at the same point where the peak of the waveform is located. For each received waveform, we first located the peak position (Tmax) of the surface return, and selected the waveform corresponding to the time interval T0 to T1 (solid red line in Fig. 4) to realize the exponential fitting; then, the exponential function was extended inversely (dashed red line in Fig. 4) to the position of Tmax, where the function was considered to have the value of Pc. Finally, Pc was subtracted from Pmax to obtain Ps; thus, Ps/Pc was obtained.
Fig. 4. Schematic of calculation of specular reflection energy (Ps) and volume backscatter energy (Pc). The attenuation process is represented as the waveform from T0 to T1. Pc is obtained from an exponential function at the peak position and Ps is obtained by subtracting Pc from Pmax.
For very high Ps/Pc, the specular reflection dominates the surface return, whereas for small Ps/Pc, the specular reflection is weaker than the volume backscatter and the peak energy is mainly contributed by the volume backscatter. Ps/Pc below 0.0 does not imply a “negative” specular reflection but that the peak of the waveform is almost entirely contributed by Pc and Ps is so weak that it is possibly submerged at the leading edge of the surface return waveform. Different Ps/Pc values affect the full waveform processing algorithms, where the peak position will alternate randomly back and forth at a point between the specular reflection and volume backscatter [3].
3. Water surface determination algorithm
3.1 Leading edge detection algorithm
The leading edge detection algorithm is a low-complexity method for water surface detection. It has been demonstrated that the leading edge detection method has significant capability to mitigate the multipath effect but at the expense of noise performance [8][17]. The key of this method is estimating the position of zero amplitude (P0). We proposed a new method based on a simple LLE, in which two points P1 and P2 on the leading edge are selected to obtain the point of zero amplitude through linear extrapolation. As shown in Fig. 5, the parameter α determines the threshold amplitude αA, where A is the maximum amplitude of the waveform. The intersection between the threshold amplitude and waveform is the position of the first point (P1). The typical value of α for sufficiently eliminating noise at the first point and above the curvature distortion is 0.15. The second point (P2) is determined by another parameter τ, which is directly related to the time interval (tint) between the maximum amplitude (tmax in Fig. 5) and the first point of the waveform (P1 in Fig. 5). The typical value of τ is within 0.25–0.50 of tint. In practice, τ is chosen as a fixed value of 0.25. This separation reduces the effect of multipath but is highly sensitive to noise. Finally, the position of zero amplitude is obtained from the intersection between the line P1P2 and the x-axis.
Fig. 5. (a) Linear-approximation-based leading edge algorithm (LLE). (b) Waveforms from Mapper5000 and CZMIL. The LLE results are shown as red symbols.
3.2 Other algorithms for comparison
3.2.1 Peak and HPP position algorithm
The peak position detection method (PD) focuses on searching for the deflection point of the waveforms but makes no special consideration of the radiometric features [18]. We adopted a one-dimensional signal processing method proposed by Wong and Antoniou to detect the peak position [19]. This method uses a finite-duration impulse response (FIR) zero-phase digital filter to suppress the noise and preserve the desired amount of information. The FIR zero-phase digital filter is expressed as follows:
3.2.2 SVB algorithm
Another type of surface extraction method is the mathematical approximation, which uses a combination of mathematical formulas to fit the received waveform. These formulas usually consider the waveform characteristics and laser radiative transfer process. Several fitting methods have been developed in the previous studies to approximate the water column response [21–24]. Schwarz et al. [6][9] proposed the SVB algorithm, wherein the received waveform pm(t) can be defined as the convolution of the system waveform hm(t) and the differential backscatter cross section σm(t). pm(t) is expressed as
where σm(t) for a water body is represented as two boxcar-shaped functions, one for the surface and one for the bottom, in addition to one exponential function segment for the water column (Fig. 6(a)). The SVB model generally considers that the water column is homogeneous and its exponential function segment depends on the parameter γ. Then, 10 unknown parameters E0…E3, τ0…τ4, and γ (Fig. 6(a)) related to the composite σm(t) were solved using the parameter optimization procedure described by Schwarz et al. [9]. From the set of parameters, τ0 was used to obtain the range to the water surface.Fig. 6. (a) SVB model with 10 parameters E0…E3, τ0…τ4, and γ. Surface and bottom are modeled as boxcars to account for pulse broadening. Effects from the water volume are modeled as an exponential segment. (b) The Mapper5000 waveform shows excellent agreement with the SVB fitting waveform.
More importantly, the system waveform hm(t) of Mapper5000 is very different from that of RIEGL VQ-880-G used by Schwarz et al. [6][9]. However, it can be defined as an exponentially modified Gaussian function (EMG) [25], which is similarly obtained via the convolution of a standard Gaussian function and an exponential function as
4. Results
4.1 Apparent water surface location uncertainty
Different contributions of the specular reflection and volume backscatter return have significant effects on the AWSL. Figure 7(a) and 7(c) show the AWSLs of East Island and Dazhou Island, extracted by the PD method, in which the point color denotes the value of Ps/Pc. Two features can be clearly observed in both experimental sites. One is the two-layer structure of the point cloud, and the other is the different distribution of the points with large Ps/Pc values. The instantaneous WSH cannot be divided into two layers, but the observed height difference between the two layers is approximately 0.5 m, which can reach 1.3 m in some extreme cases. The reason for these unexpected characteristics is the shift in the peak position of the surface return, caused by the difference in the contribution from the specular reflection returns. It can be found that the top layer mainly consists of returns with strong specular reflection whereas those at the bottom layer are dominated by volume backscatter. Furthermore, the returns dominated by specular reflection (high Ps/Pc) are concentrated on the two sides of the track in Fig. 7(a), whereas they are assembled in the middle in Fig. 7(c). These differences are associated with the scanning trajectories of Mapper5000, which have a pseudo-elliptical shape due to the periodic change in the incident angle from approximately 10.5° to 15°. The installation direction of Mapper5000 was rotated 90° horizontally in these two flight experiments. The different scanning trajectories of East Island and Dazhou Island are shown in Fig. 7(b) and 7(d), respectively. The largest incident angle of 15° corresponds to the vertexes of the major axis (red circles), whereas the vertexes of the minor axis (green circles) contain the smallest incident angle of 10.5°. The high Ps/Pc points are more likely to have occurred in the region with small incident angles, which is consistent with the Cox–Munk theory [5] that a smaller incident angle has greater possibility to generate specular reflection than a larger incident angle when the water surface is rough. Thus, volume backscatter returns superimposed with different specular reflections may cause significant uncertainty in the AWSL.
Fig. 7. (a), (c) Water surface extraction results of East Island and Dazhou Island, respectively. Ps/Pc is denoted in different colors. (b), (d) Schematics of scanning trajectories corresponding to two flight missions.
Figure 8 shows the SVB decomposition results for the waveforms with different Ps/Pc values, which can more clearly illustrate the peak position shift with Ps/Pc. The signals of volume backscatter return are basically the same in Fig. 8, with the maximum magnitude of approximately 120. However, there is an obvious deviation (nearly 4.0 ns) between the peak positions of specular reflection and volume backscatter return. In the top row of Fig. 8, the surface return is almost dominated by the specular reflection, causing the AWSL to almost coincide with the peak position of specular reflection; further, the AWSL is probably much closer to the true water surface location. As Ps/Pc decreases below 2.0 (Fig. 8(d)-(f)), the AWSL gradually deviates from the peak position of specular reflection and the time difference between them can exceed 1.0 ns. When the surface return is mainly contributed by the volume backscatter, as seen in Fig. 8(g)–8(i), the AWSL shows a greater backward shift from the peak position of specular reflection, and the time bias can exceed 4.0 ns, meaning the range bias can be more than 0.6 m, which is unacceptable. Therefore, the AWSL and peak position of surface return cannot precisely indicate the actual WSH. The situation is more serious in the case of dominant volume backscatter return, where the retrieved range in air is generally evidently overestimated and hence the WSH is underestimated. The influence of the difference in the contributions of specular reflection and volume backscatter on the WSH determination should be eliminated.
Fig. 8. SVB decomposition results of waveforms with decreasing Ps/Pc. The value of Ps/Pc is (a)-(c) greater than 2, (d)-(f) between 0.5 and 2.0, and (g)-(i) less than 0.5. The three rows correspond to three types of surface returns dominated by specular reflection, volume backscatter, and both.
4.2 Effect of water optical properties
Variations in the water optical properties definitely affect the AWSL, particularly when there is no overlap between the specular reflection return and surface return. The typical volume-backscatter-dominated surface return waveforms with low Ps/Pc were selected from the three test regions, as shown in Fig. 9(a); the normalized waveforms are shown in Fig. 9(b). The three coastal regions had very different optical properties of water; hence, their LiDAR attenuation coefficient (KLiDAR) varied from 0.08 m−1 to 0.29 m−1. Although there was no synchronous field IOP data in these test areas, it can be roughly estimated that c ranged from 0.18 m−1 to 1.25 m−1 according to the approximate relationship between Kd and c proposed in [26]. The higher the value of KLiDAR is, the stronger the magnitude of surface return will be (Fig. 9(a)); more importantly, the slope of the leading edge becomes shaper with higher KLiDAR, and the peak position shifts obviously to the right with decrease in KLiDAR (Fig. 9(b)) owing to the exponential decay behavior of the volume backscatter. The time difference between the peak positions for Ganquan and East Island could reach 3.9 ns; thus, it is considered that the clearer the water is, the greater will be the deviation of the AWSL from the true water surface location.
Fig. 9. Measured surface return waveforms dominated by the volume backscatter in three flight experiments covering different IOP waters. (a) Original waveforms. (b) Normalized waveforms. Note that the peak position differences (in nanoseconds) are presented.
To further quantify the effect of water optical properties, a controllable tank experiment (see section 2.2) was conducted and the uncertainties of the four WSH determination algorithms (PD, HPP, SVB, and LLE) were fully evaluated in waters with a wide range of IOPs. The return waveform from the plastic plate, which is very similar to hm(t) and can be considered as the reference water surface return to indicate the exact WSH, is shown in Fig. 10(a). The extraction results obtained with the plate waveform by using the four algorithms are also presented, which can be considered as the reference range in air for the corresponding algorithms (Lref). The range bias for each algorithm was obtained by subtracting its Lref from the derived range in air (L), which can be expressed as
Fig. 10. (a) Waveform of plastic plate and water surface location extracted from the three methods. (b), (c) Range biases of the four methods for different values of c and KLiDAR.
As seen in Fig. 10(b), the values of c of the tank water during the experiment ranged from 0.36 to approximately 6.0 m−1, which covered most of the near shore/offshore areas suitable for ALB. All the four algorithms overestimated the range in air and yielded positive biases. These biases showed an exponential downward trend with increasing c. The behaviors of the range bias in less turbid and more turbid water suggested that a significant amount of water volume backscattering in water with greater turbidity skewed the surface returns toward the actual water surface. It can be clearly found that the PD algorithm yielded the highest range bias, which could be more than 0.5 m when c < 1.0 m−1. The bias of the HPP method was the second highest, i.e., more than 0.3 m when c < 1.0 m−1. Noticeably, the LLE and SVB algorithms yielded almost the same results and had a much smaller range bias. Their maximum bias did not exceed 0.13 m. As c increased, the range bias decreased gradually and theoretically tended to be zero. However, the PD method still had a range bias of approximately 0.25 m when c was up to 6 m−1, and the range bias of the HPP method was near 0.2 m at the same point. The LLE and SVB methods had the minimum bias, which converged to approximately 0.06 m; this value can be ignored in practical bathymetry measurements. The relationship between KLiDAR and the range bias is also presented in Fig. 10(c). In the case of water with KLiDAR of near 0.08 m−1, the maximum depth penetration capability could reach 50 m and the range bias in air was just below 0.1 m. The remained bias may be induced by the magnitude of volume backscatter returns and signal-to-noise level of the LiDAR detector.
These results confirmed the suggestion by Schwarz et al. [6] that the actual location of the water surface could be determined from the leading edge of the observed waveforms. Thus, similar to the SVB, the LLE method can effectively suppress the variation in the volume backscatter returns and greatly reduce the range bias; moreover, it is comparatively much simpler, more robust, and requires much less computation time, as discussed in the following section.
4.3 Assessment with different rough surfaces
The LLE method, which has been proven the least susceptible to variation in the IOPs, should be further evaluated using waveforms with different Ps/Pc values based on different flight experiment datasets from Mapper5000; the performances of the PD, HPP, and SVB algorithms were similarly assessed. The WSH results of the four methods for East Island are presented in Fig. 11, in which the results of the NIR channel (HNIR) are presented as red dots. Because the NIR channel waveform contains only the specular reflection from the water surface, the HNIR results of the PD algorithm are considered as the reference WSH. For each algorithm, its WSH (denoted as Hm was derived from L as follows:
where θ is the incident angle of the output laser beam, and devm refers to the systematic displacement between LNIR and Lm (Fig. 10(a)), which is calculated from only the specular-reflection dominated scanning points with very high Ps/Pc.Fig. 11. Hm derived from four algorithms in East Island: (a) Peak detection, (b) Half peak power, (c) SVB, and (d) Leading edge detection. The color of the cloud points indicates the value of Ps/Pc; the HNIR results are presented in (a) as red points. The pink circles marked A and N refer to the scanning points with abnormal and normal profiles, respectively.
Except for the PD method, the other three methods show great improvement in reducing the water surface continuity. However, in Fig. 11(b), the WSH derived from the volume-backscatter-dominated waveforms (low Ps/Pc) by using the HPP algorithm appears to remain lower than that derived from the specular-reflection-dominated waveforms (high Ps/Pc). The SVB algorithm (Fig. 11(c)) appeared to successfully eliminate the influence of the different levels of contribution of specular reflection and volume backscatter, but it still depicted some dislocation points with low Ps/Pc located at the lower level. This was probably because the volume-backscatter-dominated waveforms theoretically did not meet the original setting of the SVB algorithm (Fig. 8) that all surface return waveforms should consist of inter-surface and volume. Comparatively, the LLE method (Fig. 11(d)) yielded the most uniform water surface. Figure 12(a)-(h) presents the side views of the water surface point clouds along and across the track near the East Island. In both cases, the WSH from the LLE method was the most consistent with the NIR results; additionally, for the points with dominant inter-surface reflection, which were mainly distributed at both ends of the cross-track (Fig. 12(e)-(h)), the results of all methods showed great consistency with the NIR results. According to the histograms of the WSH distributions, where the heights were subtracted from the mean of HNIR (Fig. 12(i)-(l)), the HNIR results near the East Island approached a Gaussian distribution from −0.25 to 0.25 m, indicating that the sea condition near the East Island was very calm during the flight experiment. The good consistency between the LLE and NIR results can be more clearly observed in Fig. 12 l, whereas the PD and HPP results have many points on the left of the NIR distribution, indicating that a considerable part of their WSHs were obviously underestimated. More detailed statistics for all scanning points are presented in Table 3, where the deviation (err.) between the mean WSH obtained with the LLE method and the NIR result is only −0.08 m; this value is much lower than −0.43 m for PD, −0.20 m for HPP, and −0.13 m for SVB. Moreover, the standard deviation (Std.) for the WSH obtained with LLE (0.140 m) was closest to the NIR result (0.108 m).
Fig. 12. Side views of water surface point cloud along (a-d) and across (e-h) the track near East Island. The WSH result from the NIR channel (HNIR) is also presented. (i-l) Histograms of HNIR and H derived using the various algorithms. Note that the HNIR distributions are the same in all subfigures.
Table 3. Statistics of WSH results in the three test areas
Another point of comparison was computational time. To enable meaningful run-time comparisons, all tests were performed on the same computer: a Lenovo workgroup, with two Intel Xeon CPU E5-2680 v4 (2.40 GHz processors), 64.0 GB RAM, running Windows 7. Noticeably, the comparison of the results in Table 3 also shows that LLE can greatly improve the computational efficiency. Although the SVB method also obtains reliable WSH results, its computational time is almost 100 times that of the LLE method.
Similar results were observed near Dazhou Island and Ganquan Island. For brevity, only the H results of Dazhou Island are shown in Fig. 13. Unlike in the case of East Island, the observations near Dazhou Island were affected by the waves. Owing to the rougher surface, a greater number of scanned points contained inter-surface reflection, thus having higher Ps/Pc values. The HNIR result near Dazhou Island also approached a Gaussian distribution, but had a larger variation from −0.6 to 0.7 m (Fig. 14) and larger Std. of 0.24 m (Table 3) whereas the corresponding value for East Island was 0.10 m. Although the WSH was still underestimated, the mean WSHs extracted from the green channel by the four methods had smaller deviations from the NIR result. The LLE method continued to yield the best performance; its difference from the mean WSH of NIR was only 0.01 m and its Std. was 0.27 m, which was closest to 0.24 m of the NIR method; hence, its distribution was almost identical to that of the NIR method. It is reasonable that the greater contribution of specular reflectance to the surface return results in more accurate determination of the WSH by all algorithms, which can also be proved based on the observations at Ganquan Island (Fig. 15) according to the statistical results in Table 3. Obviously, the sea surface during the experiment near Ganquan Island was rougher (Fig. 15), as the Std. of the NIR WSH was up to 0.37 m, and the mean WSHs and Std. from the four methods were almost identical to those of the NIR method (Fig. 16). Despite this, we can conclude that our LLE algorithm is a low-complexity practical method and shows the best performance in the determination of the WSH from the single green channel in situations where either specular reflection or volume backscatter return is dominant.
Fig. 13. H derived from four algorithms in Dazhou Island. (a) Peak detection. (b) Half peak power. (c) SVB. (d) LLE. The color of the cloud points indicates the value of Ps/Pc and the HNIR results are presented in Fig. 13(a) as red points. The pink circles marked A and N refer to the scanning points with abnormal and normal profiles, respectively.
Fig. 14. Histogram of HNIR and H derived using different algorithms in the Dazhou Island. (a) Peak detection, (b) Half peak power, (c) SVB, and (d) LLE. Note that the HNIR distributions are the same in all subfigures.
Fig. 15. H derived from four algorithms in Dazhou Island. (a) Peak detection. (b) Half peak power. (c) SVB. (d) LLE. The color of the cloud points indicates the value of Ps/Pc and the HNIR results are presented in Fig. 15(a) as red points.
Fig. 16. Histogram of HNIR and H derived using different algorithms in the Ganquan Island. (a) Peak detection, (b) Half peak power, (c) SVB, and (d) LLE. Note that the HNIR distributions are the same in all subfigures.
In addition, SVB is a good method to determine the water surface from the volume echo alone, but it can give unexpected results when dealing with some “abnormal” waveforms such as those from East Island (Fig. 11(c)) and Dazhou Island (Fig. 13(c)). The reason why the waveform is considered abnormal is that their waveform shape is very different from the adjacent normal waveform (Fig. 17(a)-(c)), and the surface positions extracted from these waveforms by using the PD, HPP, and SVB methods are very uncertain; in contrast, the surface positions extracted with LLE are reliable. Further, the mechanism behind these abnormal waveforms is not clear; they might be associated with turbulent flows and fluctuations near the water surface characterized by a wide range of coexisting scales of motion, from millimeters to meters [27]. As seen in Fig. 17(a), the A1 waveform from East Island had an unexpected flat shape, and its leading edge and the subsequent volume-backscattered signals coincided with the adjacent waveform, which had a normal surface return. There was no distinct peak in the surface return waveform of A1 in which the specular reflection was almost non-existent, resulting in failure of the SVB method (Fig. 17(d)). In the case of the abnormal waveform of A2 (Fig. 17(e)), although the convoluted waveform from SVB and the measured result were consistent near the peak position, the surface position remained incorrect; particularly, volume backscattering does not obey the exponential decrease and depicts an obvious increase near 70 ns. This probably means that the associated optical properties of water are vertically inhomogeneous. Similar waveforms and SVB results can be found for Dazhou Island in Fig. 17(f). These results can be attributed to the assumptions of the SVB method that the characteristics of the water surface and water body are invalid in some cases. Comparatively, the LLE method does not require any presuppositions; hence, it can be considered more universal and of low complexity.
Fig. 17. (a-c) “Abnormal” surface return waveforms and their neighboring normal waveforms from East Island and Dazhou Island (their positions are marked in Figs. 11 and 13). The water surface positions detected using the PD, HPP, and LLE methods are shown. (d-f) SVB model results and convolution of SVB model with system waveform.
4.4 Performance in shallow water
The LLE algorithm can provide a continuous and unified water surface for waters ranging from deep to shallow depths (Fig. 18(a)). In the region of very shallow depth near the coastal lines, the overlap between the bottom and surface returns may be a problem in water surface determination. The FWHM of the laser pulse width was only 0.95 to 1.2 ns, but owing to the limited bandwidth of PMT (∼400 MHz), the actual FWHM of the system waveform recorded by our LiDAR was more than 5 ns. At depths of less than 1 m, the overlap between the bottom and surface returns generally causes the peak of the surface return to be indistinct and even disappear; then, the peak of the bottom return will be mistakenly recognized as belonging to the surface return. In such cases, the SVB method has been proven successful in discriminating between returns from the surface and bottom, but the PD and HPP algorithms are unable to provide correct results. Although our LLE method also relies on the peak position of the surface return to determine P2, the τ value of 0.25 tint appears to be effective for very shallow waters (Fig. 18(b)-(d)). This is probably because the τ value of 0.25 tint in such cases may be equivalent to 0.5 tint in the case without overlap, which is still within the threshold suitable for the LLE method.
Fig. 18. (a) Example profile showing the water surface and bottom points in very shallow water. (b-d) Waveforms at depths of 0.25 m, 0.5 m, and 1.0 m. The three cases exhibit severe pulse overlapping.
5. Conclusion
Determining the WSH by using only the green laser returns is challenging. The uncertainty is mainly attributed to two causes: (1) the superposition of specular reflection and volume backscatter returns, and (2) variation in the optical properties of water. Peak position and leading edge are the most common factors for WSH determination. Two peak-position-based algorithms of PD and HPP as well as the leading-edge-based algorithm of SVB were comprehensively evaluated on the airborne Mapper5000 datasets with various surface roughness conditions and through a controllable tank experiment with a wide range of optical properties of water. The evaluation results indicated that the leading edge detection method of SVB had significantly better performance than the PD method. Uncertainties in the WSH derived using all algorithms were more exaggerated when the water volume backscattering dominated the surface return. Further, the uncertainties were sensitive to variations in the optical properties of water and exponentially increased with decreasing KLiDAR. A low-complexity leading edge method based on linear approximation (LLE) was proposed. Similar to the SVB method, the LLE method showed considerable robustness to the effect of volume backscatter returns and had a similarly small bias. The average and maximum errors of the LLE method were less than 0.06 and 0.13 m, respectively. Noticeably, the LLE has high computational speed, which is almost 100 times faster than that of the SVB method. Thus, the LLE method can be used as a standard method for extracting sea surface position from green laser (532 nm) echoes alone
Appendix: Notation
ALB | Airborne LiDAR bathymetry |
AWSL | Apparent Water Surface Location |
FOV | Field-Of-View |
FWHM | Full Width at Half Maxima |
FIR | Finite-duration Impulse Response |
HPP | Half Peak Power algorithm |
H | Water surface height derived from green channel |
HNIR | Water surface height derived from NIR channel |
IOPs | Inherent Optical Properties |
Kd | Diffuse attenuation coefficient |
KLiDAR | LiDAR attenuation coefficient |
L | Range in air |
LLE | Linear approximation of the leading edge |
Lref | The reference range in air for the corresponding algorithms |
PD | Peak Detection algorithm |
Ps | Specular reflection power |
Pc | volume backscattering power |
SVB | Surface-Volume-Bottom algorithm |
SIOM | Shanghai Institute of Optics and Fine Mechanics |
WSH | Water surface height |
Funding
Global Change and Air-Sea Interaction II Program (GASI-01-DLYG-EPAC02, GASI-01-DLYG-WIND02); Southern Marine Science and Engineering Guangdong Laboratory (Guangzhou) (GML2019ZD0602); National Natural Science Foundation of China (42276200); National Major Scientific Equipment and Equipment Development Special Task (2013YQ120343); State Key Laboratory of Satellite Ocean Environment Dynamics, Second Institute of Oceanography, MNR (SOEDZZ2103); Program of Innovation 2030 on Smart Ocean, Zhejiang University (129000*194232201, 129000+194432201).
Acknowledgments
We thank Yan He and Kaipeng Li from SIOM for providing the airborne LiDAR Mapper5000 and shipborne LiDAR. We also thank Kuifeng Luan and Zhenge Qiu from Shanghai Ocean University for helping us to conduct the tank experiment.
Disclosures
The authors declare no conflicts of interest.
Data availability
The data underlying the results presented in this paper are not publicly available at this time but may be obtained from the authors upon reasonable request.
References
1. G. C. Guenther, “Airborne lidar bathymetry,” Digital elevation model technologies and applications: the DEM Users Manual 2, 253–320 (2007).
2. G. Mandlburger, M. Pfennigbauer, and N. Pfeifer, “Analyzing near water surface penetration in laser bathymetry—A case study at the River Pielach,” ISPRS Ann. Photogramm. Remote Sens. Spatial Inf. Sci. II-5/W2, 175–180 (2013). [CrossRef]
3. G. C. Guenther, “Wind and nadir angle effects on airborne lidar water “surface” returns,” in Proc. SPIE (Ocean Optics, 1986), 637, pp. 277–286.
4. G. C. Guenther, P. E. LaRocque, and W. J. Lillycrop, “Multiple surface channels in scanning hydrographic operational airborne LiDAR Survey (SHOALS) airborne LiDAR,” in Proc. SPIE (International Society for Optics and Photonics, 1994), 2258, pp. 422–430.
5. C. Cox and W. Munk, “Measurement of the roughness of the sea surface from photographs of the sun’s glitter,” J. Opt. Soc. Am. 44(11), 838–850 (1954). [CrossRef]
6. R. Schwarz, N. Pfeifer, M. Pfennigbauer, and A. Ullrich, “Exponential decomposition with implicit deconvolution of lidar backscatter from the water column,” PFG 85(3), 159–167 (2017). [CrossRef]
7. Z. Pan, C. Glennie, P. Hartzell, J. C. Fernandez-Diaz, C. Legleiter, and B. Overstreet, “Performance assessment of high resolution airborne full waveform LiDAR for shallow river bathymetry,” Remote Sens. 7(5), 5133–5159 (2015). [CrossRef]
8. V. Dizdarevic and K. Witrisal, “Statistical UWB range error model for the threshold leading edge detector,” in 6th International Conference on Information, Communications & Signal Processing (IEEE, 2007), pp. 1–5.
9. R. Schwarz, G. Mandlburger, M. Pfennigbauer, and N. Pfeifer, “Design and evaluation of a full-wave surface and bottom-detection algorithm for LiDAR bathymetry of very shallow waters,” ISPRS Journal of Photogrammetry and Remote Sensing 150, 1–10 (2019). [CrossRef]
10. W. Chen, Y. Lu, J. Li, H. Jian, X. Cui, S. Yin, Y. Wang, and C. Chu, “Developing of laser airborne depth mapping in China,” Proc. Third International Asia-Pacific Environmental Remote Sensing Remote Sensing of the Atmosphere, Ocean, Environment, and Space (SPIE, 2003), 4893, pp. 414–421.
11. Y. He, S. Hu, W. Chen, X. Zhu, Y. Wang, Z. Yang, X. Zhu, D. Lü, J. Yu, T. Huang, X. Xi, S. Qu, and B. Yao, “Research progress of domestic airborne dual-frequency LiDAR detection technology,” Laser Optoelectron. Prog. 55(8), 082801 (2018). [CrossRef]
12. G. Tuell, K. Barbor, and J. Wozencraft, “Overview of the coastal zone mapping and imaging lidar (CZMIL): a new multisensor airborne mapping system for the US Army Corps of Engineers,” Proc. SPIE 7695, 76950R (2010). [CrossRef]
13. E. Fuchs and G. Tuell, “Conceptual design of the CZMIL data acquisition system (DAS): integrating a new bathymetric lidar with a commercial spectrometer and metric camera for coastal mapping applications,” Proc. SPIE 7695, 76950U (2010). [CrossRef]
14. G. C. Guenther and R. W. Thomas, “Effects of propagation-induced pulse stretching in airborne laser hydrography,” Proc. SPIE 0489, 287–296 (1984). [CrossRef]
15. G. C. Guenther, R. W. Thomas, and P. E. LaRocque, “Design considerations for achieving high accuracy with the SHOALS bathymetric lidar system,” Proc. SPIE 2964, 54–71 (1996). [CrossRef]
16. K. Li, Y. He, P. Chen, Y. Zhang, F. Liu, X. Zhu, and W. Chen, “Miniaturized underwater polarized radiation measuring instrument,” Proc. SPIE 11191, 1119115 (2019). [CrossRef]
17. I. Sharp, K. Yu, and Y. J. Guo, “Peak and leading edge detection for time-of-arrival estimation in band-limited positioning systems,” IET Commun. 3(10), 1616–1627 (2009). [CrossRef]
18. C. Wang, Q. Li, Y. Liu, G. Wu, P. Liu, and X. Ding, “A comparison of waveform processing algorithms for single-wavelength LiDAR bathymetry,” ISPRS Journal of Photogrammetry and Remote Sensing 101, 22–35 (2015). [CrossRef]
19. H. Wong and A. Antoniou, “One-dimensional signal processing techniques for airborne laser bathymetry,” IEEE Trans. Geosci. Remote Sensing 32(1), 35–46 (1994). [CrossRef]
20. V. Ramnath, V. Feygels, Y. Kopilevich, J. Y. Park, and G. Tuell, “Predicted bathymetric lidar performance of coastal zone mapping and imaging lidar (CZMIL),” Proc. SPIE 7695, 769511 (2010). [CrossRef]
21. W. Wagner, A. Ullrich, V. Ducic, T. Melzer, and N. Studnicka, “Gaussian decomposition and calibration of a novel small-footprint full-waveform digitising airborne laser scanner,” ISPRS Journal of Photogrammetry and Remote Sensing 60(2), 100–112 (2006). [CrossRef]
22. J. Reitberger, P. Krzystek, and U. Stilla, “Analysis of full waveform lidar data for tree species classification,” International Archives of Photogrammetry, Remote Sensing and Spatial Information Sciences 36(Part 3), 228–233 (2006).
23. H. Abdallah, J.-S. Bailly, N. N. Baghdadi, N. Saint-Geours, and F. Fabre, “Potential of space-borne LiDAR sensors for global bathymetry in coastal and inland waters,” IEEE J. Sel. Top. Appl. Earth Observations Remote Sensing 6(1), 202–216 (2013). [CrossRef]
24. L. Abady, J.-S. Bailly, N. Baghdadi, Y. Pastol, and H. Abdallah, “Assessment of quadrilateral fitting of the water column contribution in lidar waveforms on bathymetry estimates,” IEEE Geosci. Remote Sensing Lett. 11(4), 813–817 (2014). [CrossRef]
25. H. Wong and A. Antoniou, “Characterization and decomposition of waveforms for LARSEN 500 airborne system,” IEEE Trans. Geosci. Remote Sensing 29(6), 912–921 (1991). [CrossRef]
26. J. G. Shannon, “Correlation of beam and diffuse attenuation coefficients measured in selected ocean waters,” Proc. SPIE 0064, 3–11 (1975). [CrossRef]
27. D. J. Bogucki, J. Piskozub, M.-E. Carr, and G. D. Spiers, “Monte Carlo simulation of propagation of a short light beam through turbulent oceanic flow,” Opt. Express 15(21), 13988–13996 (2007). [CrossRef]
