Shallow Water Measurements Using a Single Green Laser Corrected by Building a Near Water Surface Penetration Model

Abstract

To reduce the size and cost of an integrated infrared (IR) and green airborne LiDAR bathymetry (ALB) system, and improve the accuracy of the green ALB system, this study proposes a method to accurately determine water surface and water bottom heights using a single green laser corrected by the near water surface penetration (NWSP) model. The factors that influence the NWSP of green laser are likewise analyzed. In addition, an NWSP modeling method is proposed to determine the relationship between NWSP and the suspended sediment concentration (SSC) of the surface layer, scanning angle of a laser beam and sensor height. The water surface and water bottom height models are deduced by considering NWSP and using only green laser based on the measurement principle of the IR laser and green laser, as well as employing the relationship between NWSP and the time delay of the surface return of the green laser. Lastly, these methods and models are applied to a practical ALB measurement. Standard deviations of 3.0, 5.3, and 1.3 cm are obtained by the NWSP, water-surface height, and water-bottom height models, respectively. Several beneficial conclusions and recommendations are drawn through the experiments and discussions.

1. Introduction

Airborne LiDAR bathymetry (ALB) is an accurate, cost-effective, and rapid technique for shallow water measurements [1]. In general, integrated infrared (IR) and green ALB systems simultaneously and collinearly emit green and IR lasers. These two lasers are used for water bottom and water surface detection, respectively, because of their wavelength and attenuation [1,2,3].

A typical waveform of a green laser mainly includes the returns from the air–water interface, water volume backscatter, and water bottom [1,4]. If the returns in the waveform are accurate, then the green laser alone may be used in an ALB system, thereby making this system compact and convenient. However, Guenther et al. (2000) rejected such a possibility due to the water surface uncertainty of green lasers [1,3]. In the waveform of green lasers, the surface return is a linear superposition of the energy that is reflected from the actual air–water interface, as well as the energy backscattered from particulate materials in the water volume just under the interface [1]. Therefore, the first return cannot exactly represent the water surface but reflects a certain level of penetration into the water column [5]. This phenomenon is called the water surface uncertainty of the green laser [1]. Such penetration is referred to as near water surface penetration (NWSP) [5]. Due to the existence of the surface uncertainty, an additional IR laser is used in the integrated IR and green ALB systems to determine the actual water surface height and accurate water depth [2]. The integrated IR and green (i.e., double laser) ALB systems, such as the Optech Coastal Zone Mapping and Imaging LiDAR (CZMIL), AHAB HawkEye II and HawkEye III, improve the measurement accuracy, but add to the costs and weights of ALB systems due to the double laser. To make the systems compact and cost-effective, single green laser ALB systems no longer use the primary IR laser, but only emit and receive the green laser. These systems have been extensively tested in coastal and riverine measurements, such as Fugro LADS LADS-MK3, Optech Aquarius, USGS Experimental Advanced Airborne Research LiDAR (EAARL) and EAARL-B, and RIEGL VQ-820-G [6,7,8,9,10,11,12,13,14,15]. As they ignore the water surface uncertainly, the simplified systems conduct the measurements despite the loss of accuracy [5,16], and cannot comply with the requirements of highly accurate applications.

If NWSP is accurately estimated, then the surface location detected by the surface return of green laser can be corrected, and the accurate water bottom location and water depth can be determined by combining the bottom return of the green laser [1]. Therefore, a green ALB system can obtain accurate measurements. At present, theoretical analysis [2,17] and statistical analysis [5,18] are used to estimate NWSP. Guenther (1985) theoretically analyzed the generation mechanism of NWSP, applied LiDAR equation to estimate the time difference between the air–water interface return and water column return, and determined the relationship between time difference and NWSP. Time difference is undoubtedly related to hydrologic conditions and ALB system parameters [1]. The theoretical analysis provides the NWSP estimation and factors that influence NWSP. However, the transmission beam is assumed to be a triangular pulse characterized by multiple forward downward scattering, high-angle backscattering, and multiple forward scattering back to the surface [17]. In general, the assumption is inconsistent with the actual situation, such as the Gaussian pulse used in Optech CZMIL [19]. In addition, many parameters used in the LiDAR equation should be manually given. Hence, the theoretical analysis method should be improved further. Statistical analysis provides another method to estimate NWSP. Mandlburger et al. (2013) used the water surface height derived from the IR laser as a reference to statistically analyze the distributions of NWSP in ponds and rivers. The results show that NWSP ranges from 10 to 25 cm, and is high in turbid water. Moreover, it is possible to approximate the reference water surface using green surface returns by statistical distribution correction [5]. The statistical analysis method is simple and efficient in calm and clear rivers and lakes, and can accurately reflect the spatial distribution of NWSP. However, this method may be affected by the undulating water surface and turbid coastal waters.

The aforementioned contributions promote our understanding of NWSP and its influencing factors. However, the precise relationship between NWSP and its influencing factors has yet to be established. Therefore, the current study proposes a method to develop the relationship model and conduct the shallow water measurements using a corrected single green laser. The structure of this paper is as follows. Section 1 introduces the existing ALB systems and studies on NWSP. Section 2 provides the detailed method of building the NWSP model. Section 3 deduces the mathematical height models of the water surface and water bottom points. Section 4 presents the validation and analysis of the proposed method through experiments. Section 5 provides the corresponding discussion. Lastly, Section 6 presents several beneficial conclusions and recommendations that are drawn from the experiments and discussions.

2. Building the NWSP Model

2.1. Comprehensive NWSP Model

Guenther et al. (1985) used the LiDAR equation to estimate the time difference between the air–water interface and the water column returns. The time difference is associated with NWSP, as well as relating to the water turbidity, transmitted pulse width, pulse detection algorithm, beam scanning angle (i.e., angle of the incident beam to the local vertical), and laser spot diameter on the water surface [1]. Among these factors, the transmission pulse width of an ALB system is fixed, and the performance of the pulse-detection algorithm adopted in the waveform detection is known. The two effects can be regarded as constants. The hydrological parameters of seawater, such as turbidity; and the measurement parameters of the ALB systems, such as the beam scanning angle and spot diameter, mainly influence NWSP. Therefore, NWSP (i.e., Δd) can be expressed as follows:

$$\Delta d={\beta }_1\phi +{\beta }_2{\phi }^2+{\beta }_{3}F+{\beta }_4{F}^2+{\beta }_{5}Turb+{\beta }_{6}Tur{b}^2+{\beta }_7$$

(1)

where β₁–β₇ are the model coefficients; φ and F are the scanning angle and spot diameter of the laser, respectively; and Turb is the turbidity of the measured sea water.

The correlations between the turbidity and suspended sediment concentration (SSC) have been investigated through extensive experiments [20,21,22,23,24,25,26]. Although turbidity depends on SSC, as well as particle composition and size distribution, many experiments have shown that a good linear relationship exists between turbidity and SSC [20,22,24,25,26]. Therefore, Turb is expressed as follows:

$$Turb=aC+b$$

(2)

where C is SSC, and a and b are the coefficients that vary with different regions and time. However, these coefficients can be regarded as constants in the same region and in a short period [20,22,24,25,26].

F can be calculated using the sensor height H and LiDAR divergence angle γ [27] using the following equation:

$$F=H\tan \left(\gamma \right)$$

(3)

The ALB system determines the transmitted pulse characteristics (i.e., initial radius and divergence angle) [28]. Thus, γ should no longer be included in the model as an independent variable. We obtain the following equation when Equations (2) and (3) are substituted into Equation (1):

$$\Delta d={\beta }_1\phi +{\beta }_2{\phi }^2+{\beta }_{3}H+{\beta }_4{H}^2+{\beta }_{5}C+{\beta }_6{C}^2+{\beta }_7$$

(4)

Equation (4) is also the NWSP model proposed in this study.

By employing the NWSP (i.e., Δd₀), which is obtained by subtracting the known water surface height from that of the green laser, as a reference, the deviation (i.e., ε_Δd) of the NWSP model can be calculated using the following equation:

$$\begin{array}{l}{\epsilon }_{\Delta d}=\Delta d-\Delta d_0\\ \Delta d_0=h_s^g-h_s^0\end{array}$$

(5)

where h^g_s is the water surface height derived from the pulse waveforms of the green laser. h⁰_s is the known water surface height, which can be obtained using an IR laser scanner mounted on the same platform as the green laser scanner; this setup is called the Mandlburger’s calibration scheme [5].

2.2. Development of the NWSP Model

The NWSP model presented in Equation (4) can be developed using the following ALB measurement and hydrological data:

(1) Measurement data of the ALB system, such as the water surface height derived from the green laser h^g_s and that derived from IR laser h^r_s, beam scanning angle φ, and sensor height H.

(2) Hydrological data of the measured water area, such as the SSC of the surface layer C.

In the preceding data, h^g_s, h^r_s, φ, and H can be extracted from the ALB records; and C can be obtained through the on-site sampling and laboratory analysis. After obtaining these data, the NWSP equation presented in Equation (4) can be established in each of the water surface points. The matrix form can be expressed as follows:

$${\mathrm{V}}_{\mathrm{n}\times 1}={\mathrm{B}}_{\mathrm{n}\times 7}{\mathrm{X}}_{7\times 1}-{\mathrm{l}}_{\mathrm{n}\times 1}$$

(6)

where n is the point number.

$$\begin{gathered} \mathrm{B}={\left(\begin{array}{ccccccc}{\phi }_1& {\phi }_1^2& H_1& H_1^2& C_1& C_1^2& 1\\ {\phi }_2& {\phi }_2^2& H_2& H_2^2& C_2& C_2^2& 1\\ {\phi }_3& {\phi }_3^2& H_3& H_3^2& C_3& C_3^2& 1\\ \cdots & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots \\ {\phi }_n& {\phi }_n^2& H_n& H_n^2& C_n& C_n^2& 1\end{array}\right)}_{n\times 7} \\ \mathrm{l}={\left(\begin{array}{ccccc}\Delta d_1& \Delta d_2& \Delta d_3& \cdots & \Delta d_n\end{array}\right)}_{n\times 1}^{\mathrm{T}} \\ \mathrm{X}={\left(\begin{array}{ccccccc}{\beta }_1& {\beta }_2& {\beta }_3& {\beta }_4& {\beta }_5& {\beta }_6& {\beta }_7\end{array}\right)}^{\mathrm{T}} \end{gathered}$$

X can be calculated based on the least squares principle as follows:

$$\mathrm{X}={({\mathrm{B}}^{\mathrm{T}}\mathrm{B})}^{-1}{\mathrm{B}}^{\mathrm{T}}\mathrm{l}$$

(7)

2.3. Variable Selection of the NWSP Model

Theory and experience give only a general direction as to which of a pool of candidate variables (including transformed variables) should be included in the regression model. The actual set of predictor variables used in the final regression model must be determined by analysis of the data. Determining this subset is called the variable selection problem. The goal of variable selection becomes one of parsimony: achieving a balance between simplicity (i.e., as few regressors as possible) and fit (i.e., as many regressors as needed) [29]. The variable selection of the regression model can be determined using a stepwise regression procedure.

The NWSP model shown in Equation (4) is only an initial model that considers the main influencing factors and can be optimized by stepwise regression. Moreover, t-test is adopted to conduct significance tests on the regression coefficients of the NWSP model. Detailed descriptions of stepwise regression and t-test theories are appended in Appendix A.

3. Height Models of the Green ALB System

Figure 1 shows the measurement mechanism of the ALB systems. The vector model for calculating the water bottom point can be expressed as follows:

$$\begin{array}{l}\mathrm{r}=\mathrm{o}+{\mathrm{r}}_{air}\frac{\Delta t_{air}{c}_{air}}{2}+{\mathrm{r}}_{water}\frac{\Delta t_{water}{c}_{water}}{2}\\ {\mathrm{r}}_{air}=\left(\sin \phi ,\cos \phi \right)，{\mathrm{r}}_{water}=\left(\sin \theta ,\cos \theta \right)\end{array}$$

(8)

where O is the scanner origin, r_air and r_water are the unit vectors of the laser propagating in the air and water, respectively; Δt_air and Δt_water are the corresponding round trip times; and c_air and c_water are the corresponding velocities, respectively. Moreover, φ and θ are the incidence angle in the air and refracted angle in the water, respectively.

If the IR and green lasers are used for the measurement, then Equation (8) can be rewritten as follows:

$$\mathrm{r}=\mathrm{o}+{\mathrm{r}}_{air}\frac{\left(t_1-t_0\right)c_{air}}{2}+{\mathrm{r}}_{water}\frac{\left(t_3-t_1\right)c_{water}}{2}$$

(9)

where t₀–t₃ are depicted in Figure 1.

Once t₁ is replaced by the round trip time of the green surface return t₂, the coordinates of the water surface and water bottom points can be achieved only by the green laser. In the substitution, it is the key to determining the accurate time delay Δt₁₂ of the green surface return relative to the IR surface return.

$$t_1=t_2-\Delta t_{12}$$

(10)

Δt₁₂ is associated with NWSP. Figure 1 shows that Δt₁₂ can be expressed as follows:

$$\Delta t_{12}=\frac{2\Delta d}{\cos (\phi )c_{air}}$$

(11)

Thereafter, we obtain the following equation:

$$t_1=t_2-\frac{2\Delta d}{\cos (\phi )c_{air}}$$

(12)

Equation (9) shows that the vector model m of the sea surface point can be calculated as follows:

$$\begin{array}{cc}\hfill \mathrm{m}& =\mathrm{o}+{\mathrm{r}}_{air}\frac{\left(t_2-\frac{2\Delta d}{\cos (\phi )c_{air}}-t_0\right)c_{air}}{2}\hfill \\ & =\mathrm{o}+{\mathrm{r}}_{air}\frac{\left(t_2-t_0\right)c_{air}}{2}+{\mathrm{r}}_{air}\frac{-\Delta d}{\cos \phi }\hfill \\ & =\mathrm{n}+{\mathrm{r}}_{air}\frac{-\Delta d}{\cos \phi }\hfill \end{array}$$

(13)

where n is the vector model of the water surface point obtained by the green surface return.

The vector model r of the water bottom point can be calculated as follows:

$$\begin{array}{cc}\hfill \mathrm{r}& =\mathrm{o}+{\mathrm{r}}_{air}\frac{\left(t_2-\frac{2\Delta d}{\cos (\phi )c_{air}}-t_0\right)c_{air}}{2}+{\mathrm{r}}_{water}\frac{\left(t_3-t_2+\frac{2\Delta d}{\cos (\phi )c_{air}}\right)c_{water}}{2}\hfill \\ & =\mathrm{o}+{\mathrm{r}}_{air}\frac{\left(t_2-t_0\right)c_{air}}{2}+{\mathrm{r}}_{water}\frac{\left(t_3-t_2\right)c_{water}}{2}+{\mathrm{r}}_{water}\frac{\Delta d}{\cos (\phi )c_{air}}{c}_{water}-{\mathrm{r}}_{air}\frac{\Delta d}{\cos \phi }\hfill \\ & =\mathrm{u}+\mathrm{e}\hfill \end{array}$$

(14)

where u is the vector model of the water bottom point and is obtained by the green bottom return, and e is its bias.

$$\begin{gathered} \mathrm{u}=\mathrm{o}+{\mathrm{r}}_{air}\frac{\left(t_2-t_0\right)c_{air}}{2}+{\mathrm{r}}_{water}\frac{\left(t_3-t_2\right)c_{water}}{2} \\ \mathrm{e}={\mathrm{r}}_{water}\frac{\Delta d}{\cos (\phi )c_{air}}{c}_{water}-{\mathrm{r}}_{air}\frac{\Delta d}{\cos \phi } \end{gathered}$$

By combining Equations (13) and (14), the coordinates of the water surface and water bottom points relative to the scanner origin O can be respectively presented as follows:

$$\left(S_s,h_s\right)=(\frac{1}{2}\left(t_2-t_0\right)c_{air}\sin \phi -\Delta d\tan \phi ,\frac{1}{2}\left(t_2-t_0\right)c_{air}\cos \phi -\Delta d)$$

(15)

$$\begin{array}{cc}\hfill \left(S_b,h_b\right)=& (\frac{1}{2}\left(t_2-t_0\right)c_{air}\sin \phi +\frac{1}{2}\left(t_3-t_2\right)c_{water}\sin \theta -\Delta d\left(\frac{{\sin }^2\phi -{\sin }^2\theta }{\cos \phi \sin \phi }\right),\hfill \\ & \frac{1}{2}\left(t_2-t_0\right)c_{air}\cos \phi +\frac{1}{2}\left(t_3-t_2\right)c_{water}\cos \theta -\Delta d\left(1-\frac{\sin 2\theta }{\sin 2\phi }\right))\hfill \end{array}$$

(16)

where (S_s, h_s) are the coordinates of the water surface point relative to O, and (S_b, h_b) are those of the water bottom point. The other symbols are the same as those in the previous equations.

If we consider only the height, then h_s and h_b can be obtained by combining the green height and the estimation of the NWSP model:

$$\begin{array}{l}{h}_s=h_s^g-\Delta d_{Model}\\ h_b=h_b^g-\Delta d_{Model}\left(1-\frac{\sin 2\theta }{\sin 2\phi }\right)\end{array}$$

(17)

where h^g_s and h^g_b are the heights of the water surface and water bottom points, respectively, which are derived from the green laser. Δd_model is the estimation of the NWSP model.

If we have an IR laser from the IR scanner mounted on the same platform as the green laser scanner, we can use the heights derived from the IR laser as reference. The errors of h_s and h_b (i.e., ε_s and ε_b) can be expressed as follows:

$$\begin{array}{l}{\epsilon }_s=h_s-h_s^r\\ {\epsilon }_b=h_b-h_b^{rg}\end{array}$$

(18)

where h^r_s and h^rg_b represent the heights of the water surface and water bottom points, respectively, which have been derived from the IR and green lasers.

4. Experiment and Analysis

4.1. Data Acquisition

To validate the reliability and accuracy of the proposed method, an ALB measurement was conducted in December, 2014, in a high turbidity water area (i.e., 18 km × 6 km) near Lianyungang, Jiangsu Province, China. In this measurement, the weather was sunny, the sea surface was calm, and wind speed was below 10 km/h. Accordingly, the ALB data were collected using CZMIL (see the primary technical parameters in

Table 1. Technological parameters of the CZMIL system.

Performance Index	Parameter
Operating altitude	400 m (nominal)
Aircraft speed	140 kts (nominal)
Pulse width	2.2 ns
Pulse repetition frequency	10 kHz
Circular scan rate	27 Hz
λ	IR: 1064 nm; green: 532 nm
Maximum depth single pulse	Kd·Dmax = 3.75–4.0 daytime (bottom reflectivity >15%)
Minimum depth	<0.15 m
Depth accuracy	(0.32 + (0.013 depth)²)½ m, 2σ
Sounding scope	0–30 m
Horizontal accuracy	(3.5 + 0.05 depth) m, 2σ
Sounding density	2 m × 2 m nominal
Scan angle	20° (fixed off-nadir, circular pattern)
Swath width	294 m (nominal)

). The CZMIL system emits green (λ = 532 nm) and IR (λ = 1064 nm) laser pulses in a collinear manner [19]. The beam scanning pattern is circular with a fixed angle [30]. The diffuse attenuation coefficient (Kd) derived from the moderate resolution imaging spectroradiometer (MODIS) was above 1.5 m⁻¹. Prior to the measurement, the ground control points were set to calibrate the CZMIL system. During the measurement, five SSC sampling stations were arranged in the measurement area (see in Figure 2). Seawater samples were collected at 0.5 m below the water surface of each sampling station in situ using horizontal water samplers and submitted to the laboratory for analysis. Each water sample was filtered, dried, and weighed.

Table 2. Suspended sediment concentration (SSC) of the surface layer in each sampling station.

Sampling Station	SSC of the Surface Layer (mg/L)
1	315
2	122
3	134
4	110
5	185

lists each station’s SSC. After the field measurement of 16 lines, the CZMIL data were processed using the Optech HydroFusion software, and the 3D point cloud data of the water surface and water bottom points were obtained. Figure 2 shows the data.

4.2. Construction and Optimization of the NWSP Model

By considering the SSC sampling station as the center, a representative water area of 100 m × 100 m is selected. A total of 16,076 water surface point pairs of IR and green lasers are obtained in the five representative water areas. By using the water surface height of the IR laser as reference, NWSP Δd in each point pair can be calculated using Equation (5). In addition, the sensor height H and beam scanning angle φ in each measuring point are extracted from the raw ALB records, and the SSC of the surface layer C in each sampling station is also obtained from Table 2.

Table 3. Statistical parameters of the four-type data used in the modeling.

	SSC of the Surface Layer C (mg/L)	Sensor Height H (m)	Beam Scanning Angle φ (°)	NWSP Δd (cm)
Max.	315	438	23.5	39.4
Min.	110	408	17.2	17.2
Median	134	423	20.1	28.5

lists the statistical results of the four-type data. Evidently, Δd ranges from 17.2 cm to 39.4 cm and has a mean of 28.6 cm and standard deviation of 3.8 cm. This change implies that the effect of NWSP on the water surface height derived from the green laser is significant. Thus, an NWSP model should be built to improve the water surface height of the green laser.

The four-type data of 16,076 point pairs are used for statistical analysis to assess the change of Δd with the beam scanning angle φ, sensor height H and SSC of the surface layer C. Figure 3 lists the statistical relationships. Evidently, the relationships of Δd varying with φ and H are approximately linear (see Figure 3a,c), and Δd changing with C is nonlinear (see Figure 3e). These relationships show that the proposed empirical model of Δd (see Equation (4)) is reasonable.

A total of 14,290 of the point pairs are used to build the NWSP model shown in Equation (4) by Equations (6) and (7). The rest are applied to test the model. The model coefficients β₁–β₇ are estimated via the regression analysis (see

Table 4. Coefficients and their significances to the comprehensive and optimized models. SE: standard error; t: t-statistics.

Item	Coefficient (Units)	Comprehensive Model				Optimized Model
		Value	SE	t	p	Value	SE	t	p
φ	β1 (m·deg−1)	8.17 × 10−3	6.30 × 10−3	1.2981	0.1942	8.44×10−3	2.80×10−4	30.093	0.0000
φ2	β2 (m·deg−2)	1.66 × 10−6	1.54 × 10−4	0.0107	0.9913
H	β3 (1)	8.03 × 10−3	6.10 × 10−3	1.3160	0.1882
H2	β4 (m−1)	−9.74 × 10−6	7.25 × 10−6	−1.3428	0.1793	−1.9×10−7	8.89×10−8	−2.1574	0.0309
C	β5 (m·mg−1·L)	2.11 × 10−3	4.12 × 10−5	51.279	0.0000	2.12×10−3	4.09×10−5	51.899	0.0000
C2	β6 (m·mg−2·L2)	−4.63 × 10−6	8.98 × 10−8	−51.489	0.0000	−4.65×10−6	8.80×10−8	−52.843	0.0000
Constant	β7 (m)	−1.7389	1.2734	−1.3656	0.1721	−5.4×10−2	1.73×10−2	−3.1341	0.0017

). To guarantee that these coefficients are significant in the model, t-test is adopted to conduct the hypothesis tests on these regression coefficients. The standard error (SE), t-statistic (t), and p-value (p) of each coefficient is calculated (see Table 4). Accordingly, comparing the p-values of these coefficients shows that all coefficients, except β₅ and β₆, are higher than the standard α = 0.05 cutoff, thereby indicating multicollinearity in the comprehensive model shown in Equation (4) and that the model cannot substantially represent NWSP. To solve this problem, stepwise regression is adopted to optimize the model. Equation (19) shows the optimized NWSP model. Relative to Equation (4), the optimized model drops φ² and H, thereby indicating that the impact of the beam scanning angle φ on NWSP can be described through a simple linear function. Upon optimization, the p-values of the remaining parameters in the optimized model are below α. Therefore, the parameters in Equation (19) are statistically significant and should be included in the model.

$$\Delta d={\beta }_1\phi +{\beta }_4{H}^2+{\beta }_{5}C+{\beta }_6{C}^2+{\beta }_7$$

(19)

To analyze the contribution of each parameter in the optimized NWSP model, the standardized coefficients $\tilde{\beta }$ of each parameter in the model are calculated (see

Table 5. Standardized coefficients $\tilde{\beta }$ of the different parameters.

Parameter in the NWSP Model	φ	H2	C	C2
$\tilde{\beta }$	2.74 × 10−1	−2.35 × 10−2	4.314	−4.144

). A detailed description of standardized coefficients theory is appended in Appendix A. The contribution of the surface layer SSC is the largest, followed by the beam scanning angle φ and ALB sensor height H.

To accurately analyze the effects of the preceding parameters on NWSP, Equation (19) and these coefficients in the optimized NWSP model are used to reflect the relationships of Δd varying with φ, H, and C. Figure 3 shows these relationships. Relative to the statistical relationship, the relationship between each parameter and Δd in the optimized model is considerably clear. The impact of φ on Δd is a linear positive correlation (see Figure 3b), that is, Δd increases with φ. The impact of H on Δd has a monotone decreasing trend (see Figure 3d). In a relatively small range of C, Δd increases with C, whereas Δd decreases gradually with C when C increases to a certain range (see Figure 3f). The phenomenon of Δd varying with C is consistent with reality. For an ALB system with a specific laser frequency or wavelength, C directly determines the surface density of seawater, thereby affecting NWSP based on Lambert’s law [31]. When the value of C is high, the laser penetration ability will decrease, and Δd will inevitably decrease. In terms of the CZMIL instruments, the turning point of NWSP changing with the SSC of the surface layer appears at C as 220 mg/L (see Figure 3f).

4.3. Height Calculation

After obtaining the NWSP model, the heights of the water surface and water bottom points in the measured water area can be calculated using only a green laser and the NWSP model. The process is as follows:

First, the surface layer SSC at any location of the measured area is determined by an inverse distance weighting (IDW) interpolation using SSC and the coordinates of the SSC sampling stations. IDW is extensively used to interpolate the spatial data (e.g., SSC) and considered an effective method [32,33,34]. A detailed description of IDW theory is appended in Appendix A. The interpolation adopts the assumption that the change of the surface layer SSC among these sampling stations is gradual.

Second, NWSP Δd, is calculated by Equation (19) using φ, H, and C in each of the water surface points of the green laser.

Third, the heights of the water surface and water bottom points h_s and h_b are calculated using Equation (17).

Lastly, the accuracy of the height is assessed by comparing the calculation height with the reference. In the comparison, the reference heights of the water surface and bottom points h^r_s and h^rg_b are provided by the integrated IR and green lasers.

4.4. Accuracy Analysis

4.4.1. Accuracy Analysis for the NWSP Models

Using the data of the remaining 1786 point pairs in the representative water areas, the comprehensive NWSP model shown in Equation (4) and optimized model shown in Equation (19) are evaluated. The Δd in each point pair can be obtained by considering the beam scanning angle, sensor height and surface layer SSC in the two models. The NWSP model errors can be calculated by comparing the two NWSP models with the height difference derived from the IR and green surface returns. To perform statistical analysis for the model errors, the statistical parameters and corresponding probability density function (PDF) curves of the two model errors are listed in

Table 6. Statistical parameters of the comprehensive and optimized near water surface penetration (NWSP) model errors.

NWSP Model	Max/cm	Min/cm	Mean/cm	Std./cm
Comprehensive model	10.52	−9.01	0.1	3.1
Optimized model	10.51	−8.97	0.1	3.0

and drawn in Figure 4. The statistical parameters and distributions of the two model errors are nearly the same. The statistical result shows that the two models are nearly equivalent in accurately reflecting NWSP. Hence, the optimized model simplifies the expression and improves the accuracy of the comprehensive model. The preceding analysis also demonstrates that the modeling method proposed in this study is efficient.

4.4.2. Accuracy Analysis of the Height Models

The heights of the water surface and bottom points derived from the integrated IR and green lasers are used as references. The heights obtained by the height models depicted in Equation (17) are used to compare the reference heights. For the statistical analysis of the height errors,

Table 7. Statistical parameters of the water-surface height and water-bottom height model errors.

Statistic Parameter	Max./cm	Min./cm	Mean/cm	Std./cm
Error of the water-surface height model	19.8	−20.0	1.3	5.3
Error of the water-bottom height model	5.6	−8.7	0.7	1.3

lists the statistical parameters of the model errors and Figure 5 shows the corresponding PDF curves. Evidently, the accuracy of water bottom height is higher than that of the water surface height. Approximately 92.5% of the water surface height errors are below 10 cm, whereas 82.1% of the water bottom height errors are under 2 cm. The reason is that the accuracy of the surface layer SSC possibly affects the estimation of the water surface height. In the surveyed water area, the SSC evidently changes from 110 mg/L to 315 mg/L (see Table 2). The IDW interpolation may decrease the SSC accuracy at the SSC abrupt variations. The effect of SSC on NWSP is significant; thus, an inaccurate SSC will result in an inaccurate NWSP and water surface height. Nonetheless, the method proposed in this study achieves a standard deviation of 5.3 cm in the surface height estimation and that of 1.3 cm in the water bottom height estimations, thereby proving the good performance of the proposed method in the surveyed water area. The preceding statistical results also show that the water surface height and water bottom height can be determined accurately using only the corrected green laser.

5. Discussion

The proposed method provides a good technique to obtain accurate water surface height and water bottom height using a corrected single green laser. The following main factors determine the applications and accuracies of the proposed method.

(a) Reference water surface height

Only when the actual or reference water surface height is known can the actual NWSP of the green laser be calculated, and the NWSP model and height models be established. Therefore, the reference water surface height is the basis of the proposed method. In this study’s experiments, the reference height is provided by the IR laser because the double-laser ALB system (i.e., CZMIL) was adopted in the measurement. Although the processing will not influence the research result, it may become meaningless in the actual application of the proposed method, because the aim is to use a single green laser ALB system and not a double-laser one. To address this issue, two methods are provided as follows.

(1)

Addition of an assistant IR laser scanner

If an assistant IR laser scanner is mounted on the same platform as the green laser ALB system, the assistant IR laser can provide the reference water surface height. The assistant IR laser scanner is only used prior to the measurement. This scheme was also proposed by Mandlburger et al. [5] and validated to be efficient in various water areas.

(2)

Water level of calm water

If the measured water area (e.g., lakes) is calm, then the reference water surface height can be determined through the water level. The water level can easily be determined, such as through the use of Global Positioning System (GPS) or leveling measurement.

(b) SSC

SSC is the primary influencing factor in the NWSP model. However, the hydrological parameter can only be obtained in a few limited sampling stations. An interpolation has to be performed to calculate SSC at any location surrounded by these sampling stations. To guarantee the interpolation accuracy, the following principles need to be followed.

(1)

Set the sampling stations with a certain density to ensure that the SSC from these stations can reflect the SSC variations in the water area measured.

(2)

Set the stations at the representative locations. Only a few stations are arranged in the water area with small SSC change. By contrast, numerous stations have evident SSC change.

(c) Wave effect

As shown in Figure 6, the beam scanning angle φ is the angle of the incident beam to the local vertical. The angle between the incident laser and normal of the water surface at the laser point is known as the angle of incidence. If a wave slope γ exists, which is the angle of the normal to the local vertical, then the actual incident angle of the laser beam is equal to the sum of φ and γ.

Most wave energy is typically concentrated in wind wave.

Table 8. Wind speed and the parameters of its resulting waves.

Wind Speed in km/h	Average Height in m	Average Wavelength in m	Wave Slope Range in Degrees
20	0.3	10.6	(−3.2, 3.2)
30	0.9	22.2	(−4.6, 4.6)
40	1.8	39.7	(−5.2, 5.2)
50	3.2	61.8	(−5.9, 5.9)
60	5.1	89.2	(−6.5, 6.5)
70	7.4	121.4	(−7.0, 7.0)

shows the wind speed and its resulting wave height H and wavelength L [35]. As shown in Figure 6, the approximate range of the wave slope can be estimated using the following equation:

$$\gamma \in \left(-\arctan \left(2H/L\right),\arctan \left(2H/L\right)\right)$$

(20)

The optimized NWSP model in Equation (19) is used to estimate the effect of wave on NWSP. The effect ε_d can be expressed as follows:

$$\begin{array}{cc}\hfill {\epsilon }_d& =\Delta d_1-\Delta d_2\hfill \\ & ={\beta }_1\phi -{\beta }_1(\phi +\gamma )\hfill \\ & =-{\beta }_1\gamma \hfill \end{array}$$

(21)

where Δd₁ is the NWSP estimation obtained by ignoring the wave slope, whereas Δd₂ is that obtained by considering the wave slope.

Figure 7 shows the variations of ε_d with the wave slopes. Evidently, ε_d ranges from −2.7 cm to 2.7 cm at wind speed of 20 km/h. When the wind speed reaches 30 km/h, ε_d changes from −3.9 cm to 3.9 cm. ε_d also varies with the wind speed. In our experiment, the wind speed is below 10 km/h, and the effect of the wave on the NWSP changes from −1.2 cm to 1.2 cm, which can be ignored. However, if the wind speed is above 30 km/h, then ε_d changes from −4.0 cm to 4.0 cm and the maximum absolute effect is more than the accuracy of the NWSP model. Therefore, the wave effect should be considered by adding the wave slope to the beam scanning angle.

(d) Alternative methods

Shipboard echo sounding is a traditional and accurate bathymetry method. However, this method is constrained by its high operating cost, inefficiency, and inapplicability in shallow waters [36]. By comparison, remote sensing methods provide substantially flexible, efficient, and cost-effective means over broad areas. The remote sensing of bathymetry is broadly categorized into two: active non-imaging and passive imaging methods [37,38]. The active non-imaging method (as typified by ALB) can produce accurate bathymetric information over clear waters at a depth of up to 70 m [37].

Relative to ALB, the passive imaging method provides an economical and flexible method to obtain bathymetric data from space. Dörnhöfer et al. (2016) evaluated the potential of Sentinel-2 to retrieve water depth [39]. They learned that the retrieved water depths were highly correlated with echo sounding data (r = 0.95, residual standard deviation = 0.12 m) of up to 2.5 m (Secchi disk depth: 4.2 m), even though the water depths were slightly underestimated (Root Mean Square Error (RMSE) = 0.56 m). WASI-2D is a freely available software tool for analyzing atmospherically corrected multispectral and hyperspectral imagery in both optically deep and shallow water [39,40]. In deeper water, Sentinel-2A bands were incapable of allowing a WASI-2D-based separation of macrophytes and sediment, which leading to erroneous water depths [39]. Relative to the expensive ALB data, the Sentinel-2A data can be obtained immediately and for free. However, the range and accuracy of the retrieved depth need to be improved further.

These methods are complementary. In deep waters, the traditional echo sounding method can be used where the remote sensing method may be inefficient. In shallow waters, the active remote sensing method can be used to perform accurate bathymetry when the traditional echo sounding is inefficient. The passive imaging method can be used when budget and accuracy demand are low.

(e) Limitations

Our experiment is conducted in a high turbidity area (Kd > 1.5 m⁻¹), in which SSCs of the sampling stations range from 110 mg/L to 315 mg/L. The NWSP model is adequate to fit SSC for these SSC data. However, when SSC is extremely high, the model value may be negative; thus, the model will become inefficient. In terms of the CZMIL parameters, when SSC is above 490 mg/L, the beam scanning angle is 20°, and the sensor height is 420 m; thus, the negative value will appear (see Figure 3f). In this case, the polynomial function in the NWSP model will become inefficient, and a Gaussian function may be appropriate to fit the extreme high SSC.

6. Conclusions and Suggestion

By building the NWSP model and height models, the proposed method accurately determines the water surface height and water bottom height using only a corrected green laser. Compared with the results of the green ALB system, the proposed method remarkably improves the height accuracy. Relative to the integrated IR and green ALB system, the proposed method simplifies the ALB system in size and cost. Moreover, the proposed method achieves the standard deviations of 3.0, 5.3, and 1.3 cm in the estimations of NWSP, water surface height, and water bottom height, respectively, in the experiments.

The proposed method was tested in a coastal area and further tests should be carried out in a more riverine setting. The reference water surface height is the basis of the proposed method. Therefore, an auxiliary IR laser scanner is recommended in the green laser ALB system to calibrate the system and calculate the reference water surface height. In addition, to guarantee the accuracy of the proposed method, the density and the representativeness should be considered in setting SSC sampling stations based on the measured water area.